{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11.2马尔科夫链蒙特卡罗"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "马尔科夫链定义本身比较简单，它假设某一时刻状态转移的概率只依赖于它的前一个状态。\n",
    "$$p(z^{(m+1)}|z^{(1)},...,z^{(m)})=p(z^{(m+1)}|z^{(m)})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于⼀个特定的变量，边缘概率可以根据前⼀个变量的边缘概率⽤链式乘积的⽅式表⽰出来\n",
    "$$p(z^{(m+1)})=\\sum_{z^{(m)}}p(z^{(m+1)}|z^{(m)})p(z^{(m)})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于一个转移概率为$t(z',z)$的同质马尔科夫链来说，如果:\n",
    "$$P^*(z)=\\sum_{z'}T(z',z)p^*(z')$$\n",
    "那么P^*(z)是不变的，它的一个充分非必要条件是概率转移矩阵满足细节平衡条件$$P^*(z)T(z,z')=P^*(z')T(z',z)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造一个马尔科夫链使得所求的概率分布是不变的，然后从中采样，对于$m\\rightarrow\\infty$，概率分布$p(z^{(m)})$收敛到平稳不变的概率分布$p^*(z)$,与初始概率分布$p(z^{(0)})$无关"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "#sys.path\n",
    "sys.path.append(\"/home/leilichuan/machine learning/ML notes/Pattern-Recognition-and-Machine-Learning/\")\n",
    "import matplotlib.animation as animation\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from prml.rv import Gaussian, Uniform,multivariate_gaussian\n",
    "from prml.sampling import metropolis, metropolis_hastings, rejection_sampling, sir\n",
    "\n",
    "np.random.seed(666)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "current round 35\n",
      "[ 0.62499289  0.3125065   0.06250061]\n"
     ]
    }
   ],
   "source": [
    "trans_matrix=np.array([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]])\n",
    "p_z_0=np.array([0.3,0.4,0.3])\n",
    "p_z=[p_z_0]\n",
    "for i in range(100):\n",
    "    p_z_0=p_z_0@trans_matrix\n",
    "    p_z.append(p_z_0)\n",
    "    if np.allclose(p_z[i-1],p_z[i]):\n",
    "        print(\"current round\",i+1)\n",
    "        print(p_z_0)\n",
    "        break\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "current round 33\n",
      "[ 0.62500649  0.31249407  0.06249944]\n"
     ]
    }
   ],
   "source": [
    "trans_matrix=np.array([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]])\n",
    "p_z_0=np.array([0.7,0.1,0.2])\n",
    "p_z=[p_z_0]\n",
    "for i in range(100):\n",
    "    p_z_0=p_z_0@trans_matrix\n",
    "    p_z.append(p_z_0)\n",
    "    if np.allclose(p_z[i-1],p_z[i]):\n",
    "        print(\"current round\",i+1)\n",
    "        print(p_z_0)\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一般情况下，平稳分布$P^*(z)$和某一个马尔科夫状态转移矩阵T不满足细节平衡条件\n",
    "$$p^*(z)q_k(z|z^{(\\tau)})\\neq p^*(z^{(\\tau)})q_k(z^{(\\tau)}|z)$$\n",
    "为使该条件成立引入$A_k(z^*,z^{(\\tau)})$:\n",
    "$$p^*(z)q_k(z|z^{(\\tau)})A_k(z^*,z^{(\\tau)})=p^*(z^{(\\tau)})q_k(z^{(\\tau)}|z)A_k(z^{(\\tau)},z^*)$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "基于马尔科夫链的metropolis采样过程：\n",
    "\n",
    "* 设定状态转移矩阵T,设定转移次数阈值$n_1$,需要的样本数量$n_2$\n",
    "    \n",
    "* 从任意简单概率分布采样得到初始状态矩阵$p(z^{(0)})$\n",
    "    \n",
    "* $m=0$到$m=n_1+n_2-1$\n",
    "\n",
    "     * 从条件概率分布$q_k(z^{(m+1)}|z^{(m)})$采样得到样本$z^{(m+1)}$\n",
    "        \n",
    "     * 从均匀分布采样$u-uniform(0,1)$\n",
    "     \n",
    "     *如果$$u<min(1,\\frac{\\tilde p(z^*)}{\\tilde p(z^{(\\tau)})})$$则接受转移$z^{(\\tau)}\\rightarrow z^*$\n",
    "     * 否则不接受转移，$\\tau=max(\\tau-1,0)$\n",
    "* 样本集$(z_{n_1},...z_{n_1+n_2+1})$即为平稳分布对应的样本集"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd0XNW1+PHvntGo92LLVrUx7hhsK2BM6IHQDYEQA3Fo\nwXECKY+W+iCQRvLLei8h8EhMSUggpodAKKGZAAZjWzZy77a6LFmS1evM+f0xxZKsMpJmNKOZ/Vlr\nljV3zty7Zyzte+65+54rxhiUUkqFFkugA1BKKeV7mtyVUioEaXJXSqkQpMldKaVCkCZ3pZQKQZrc\nlVIqBGlyV0qpEKTJXSmlQpAmd6WUCkERgdpwenq6yc/PD9TmlVJqXCosLDxsjMkYql3Aknt+fj4b\nNmwI1OaVUmpcEpFib9rpsIxSSoUgTe5KKRWCNLkrpVQI0uSulFIhaMjkLiLRIrJORIpEZJuI3NdP\nGxGRB0Vkr4hsFpEF/glXKaWUN7yplukAzjHGNIuIDfhIRN4wxqzt0eZC4HjX4xTgEde/SimlAmDI\n5G6ct2pqdj21uR59b9+0BPirq+1aEUkWkUnGmEqfRquU6ldjYyPV1dV0dXUFOhQ1CjabjQkTJpCY\nmDjqdXlV5y4iVqAQmAY8bIz5tE+TLKC0x/My1zJN7ipkVTa0sWZvLV+an4XFIgGLo7GxkUOHDpGV\nlUVMTAwigYtFjZwxhra2NsrLywFGneC9Su7GGDtwkogkA/8QkbnGmK3D3ZiILAeWA+Tm5g737UoF\njc1lR7j5yQ3UNHUQY7Ny8bxJAYulurqarKwsYmNjAxaDGj0RITY2lqysLCoqKkad3IdVLWOMOQKs\nBi7o81I5kNPjebZrWd/3rzTGFBhjCjIyhrx6Vqmg9ObWSq7+0ydERVjIS4vld+/sxu4I3I3mu7q6\niImJCdj2lW/FxMT4ZHjNm2qZDFePHRGJAc4DdvZp9grwNVfVzCKgQcfbVSh6e/shVjy1kVmTEnn5\n1tO48/wZ7Klu5rUtgf1116GY0OGr/0tveu6TgNUishlYD7xtjPmXiKwQkRWuNq8D+4G9wKPAt3wS\nnVJB5oXCUiYlRbPqlkWkx0dx8QmTmD4xPuC9d6X68qZaZjMwv5/lf+zxswFu9W1oSgWXLruDj/fW\ncsmJk4i2WQGwWITvfWE633p6I68UlXPF/OwAR6mUk16hqpSXNpUcoamjmzOn9z5fdMGcTGZmJvDg\nu3vptjsCFJ1SvWlyV8pL/9ldjdUiLJ6W3mu5xSLcevY0DhxuobC4PkDRhY7rr7+elJQUfv/737Nq\n1SpuuOGGXq8bYzjppJN48sknvV7nbbfdxs033+zjSIObJnelvPSf3TUszE0hMdp2zGunuRL+xpIj\nYx1WSKmvr+ff//43jz32GE899RS33norN910U682zz33HHV1dVx77bVer/fOO+/k6aefZu/evb4O\nOWhpclfKCzVNHWwtb+TMGf2X8KbGRZKfFsumEu25j0ZKSgpVVVVceeWVrF+/nrq6Os4444xebR58\n8EGWLVuGzXbsTnYg+fn5fP7zn+eRRx7xdchBS5O7Ul74cE8NwDHj7T0tyE1hY8kRnPUFaiRSU1P5\n7//+b6677jrS0tJITU3lG9/4Bh0dHQDs3buXjz/+mKuuuqrX+5599llEpN9HY2MjAFdeeSVPP/00\nDkd4nBfR5K6UF/6zu4b0+EhmTxr4qsH5uckcbu6grL5tDCMLHSUlJdTX1/O///u/pKam8vzzz/O9\n732PlStX8vDDDwPw7rvvEhcXx4knntjrvWeccQaffPKJ5/Haa6+RmprKRRdd5LnSc/HixRw6dIgt\nW7aM+WcLhIDdQ1Wp8cLuMHywu4azZkwYdA6Z+bkpAGwqPUJOqk4FMFxFRUWA8+TnAw88AMA555zD\nBx98wJtvvsntt99OYWEhs2bNwmLp3S+dNGkSkyY5p4BobGzk3HPPZe7cuTz//POeNnPmzMFqtbJu\n3bpjdg6hSJO7UkPYWt5AfWvXoEMyADMzE4i2WdhYXM9lJ04eo+j6d9+r29he0RiQbc+enMi9l84Z\n9vuKioqIjY3lrrvu6rV8+vTprFu3DoCqqirS09P7ezsALS0tXHTRRVgsFv71r3/1mm8nIiKC5ORk\nqqqqhh3beKTJXakhuMfbP3/8wEkFIMJqYV52MptKtWJmJIqKijjjjDNIS0vrtbyqqorJk507y/b2\n9gEnSGtvb+eyyy6jpaWF9957j4SEhGPaREVF0d7e7vvgg5Amd6WGsLW8kSnpcaTHRw3ZdkFuCo9/\ntJ/2LrvnKtZAGEnPOdCKioo4/fTTey1ramrivffe4777nDeAS01N7bfn3dXVxZVXXkllZSX/+c9/\nSElJ6XcbR44cITU11ffBByE9oarUEHYfamL6xHiv2s7PTabLbthW0eDnqEJLS0sL+/btY9++fb2W\nP/jgg4DzwiaAGTNmcODAgV5t7HY711xzDbt27eKdd95hoBlna2pqaG1tZfr06X74BMFHk7tSg2jv\nsnOwtoUZmd7NrT0/NxmAjcU6NDMcW7ZsweFwsG/fPu6//35Wr17ND3/4Q+677z5WrlxJcrLzez3t\ntNMoKSmhpqbG895vfvObvP766/zyl7+kpKSEtWvXsnbt2mOqYjZs2ICIsHjx4jH9bIGiwzJKDWJv\ndTMOAzMmHjt+258JCdFkp8SwqVQvZhqOoqIiEhMTee2117jxxhv5xS9+wYwZM3j++edZsmSJp91Z\nZ51Famoqb775JsuWLcMYwzPPPENbWxtf+cpXeq3z0ksv5ZVXXvE8f/PNNznzzDOPGdMPVZrclRrE\n7kNNAMzI9G5YBpwlkRsO1vkrpJBUVFTECSecwLx58ygsLBywXWRkJF/96ld55plnWLZsWa+LlAZj\nt9t58cUXPSWW4UCHZZQaxK6qJiKtFvLS4rx+z4LcZCob2qls0IuZvFVUVOR17fldd93F6tWr2b17\nt9frf/7554mJiWHp0qUjDXHc0eSu1CB2HWriuAnx2Kze/6nMy3aOD28tD0yd+XhjjGHLli3MmzfP\nq/bZ2dk88cQTVFZ6f/crYwyPP/44ERHhM1gRPp9UqRHYXdXEyVOGVzo3LcM5hLO/phmY6IeoQou3\nQys9DbcHfs011wyrfSjQnrtSA2hs76KioZ3pmd6dTHVLirWRHh/J/poWP0Wm1NA0uSs1gD3uk6le\nVsr0NDU9nn01zb4OSSmvaXJXagA7q5zJffoIkvtxE+LYf1h77ipwNLkrNYDdVU3ERVrJTokZ9nun\npsdT19JJfUunHyJTamia3JUawK5DTUzPTEBk4Gl+BzI1w1k6uf+wDs2owNDkrlQ/jDHsqmoa0Xg7\nwHGuipl9elJVBciQyV1EckRktYhsF5FtIvLdftqcJSINIvKZ63GPf8JVamwcbu6kvrVrROPtANkp\nMdisohUzKmC8qXPvBu4wxmwUkQSgUETeNsZs79PuQ2PMJb4PUamxt6vKPe3AyJJ7hOuq1v1aMaMC\nZMieuzGm0hiz0fVzE7ADyPJ3YEoF0q5Do0vuAMdlxGk5pAqYYY25i0g+MB/4tJ+XF4vIZhF5Q0TG\n350ClOphz6EmUuMivbpBx0CmZsRTUtdKt93hw8iU8o7XyV1E4oEXge8ZY/peK7wRyDXGzAP+ALw8\nwDqWi8gGEdnQcz5mpYLNwdoWpqR7P1lYf6amx9FlN5TW6wRioeD9999HRNi6dWugQ/GKV8ldRGw4\nE/vTxpiX+r5ujGk0xjS7fn4dsInIMTecNMasNMYUGGMKBrpbilLBoKS2lbzU/u/V6a3jJrgqZqp1\naEaNPW+qZQR4HNhhjPmfAdpkutohIie71lvry0CVGivtXXYqG9vJTRtlck93TSCmte4qALzpuZ8G\nLAPO6VHqeJGIrBCRFa42VwFbRaQIeBBYaowxfopZKb8qq2/FGMgbZXJPirWRFqcTiHlr27ZtXHDB\nBaSmphIXF8esWbN4+OGHAXjttdc477zzmDBhAomJiSxatIi33nqr1/t/+tOfkp6ezqeffkpBQQEx\nMTF8/vOf58CBA1RXV3P55ZcTHx/PrFmzeO+993q9Nz8/nzvvvJOf/exnZGZmEh8fz3XXXUdDw+D3\nwnU4HDzwwANMmzaNqKgopk+fzpNPPtmrzUcffcTpp59OYmIiiYmJnHTSSTz//PM++MYGN2QppDHm\nI2DQS/SMMQ8BD/kqKKUCqbi2FYDc1NGNuYPzSlVN7t659NJLmTVrFk899RRRUVHs2rXLMxXwgQMH\nuOiii7jjjjuwWq288cYbXHjhhXzwwQecdtppnnW0trayfPly7r77buLi4vjOd77DsmXLiIqK4sIL\nL+Rb3/oWv/nNb/jyl79MaWkpsbFHd+CrVq1i2rRpPProo1RWVnL33Xfz9a9/fdBE/O1vf5snn3yS\ne+65hwULFvD2229z0003kZaWxiWXXEJjYyOXXHIJS5Ys4Z577vHMXX/kyBjcY9cYE5DHwoULjVLB\n6PEP95u87//L1DS1j3pd33+hyCy4/y0fRDWw7du3+3R9drvDVDe2G4fD4dP1DqampsYAZvPmzUO2\ntdvtpqury5x//vnmxhtv9Cy/9957DWDef/99z7KHH37YAOa+++7zLNu2bZsBzOuvv+5ZlpeXZ1JS\nUkxTU5Nn2VNPPWVExPP9rl692gBmy5Ytxhhj9uzZY0TE/OUvf+kV37Jly0xBQYExxpj169cbwDQ2\nNg7n6xj0/xTYYLzIsTr9gFJ9lNS1EhdpJS0uctTrmpoRR21LJw2tXT6IzP8cDsM1j67l1F+9y9KV\na3E4xmZ0NTU1lZycHFasWMGzzz5LdXV1r9fLysq4/vrrycrKIiIiApvNxltvvXXMrfYiIyM5/fTT\nPc+nTZsGwDnnnHPMsvLy8l7vPe+884iPP3qv3CuuuAJjDOvXr+835nfffReLxcIVV1xBd3e353Hu\nuefy2WefYbfbOe6444iPj+faa6/ln//859j02F00uSvVR3FtC7lpcSOaMKyvqa6TqvvGyUnV2pZO\nCovr6XYYCovrqR2jWS0tFgtvvfUWmZmZ3HTTTWRmZnL66aezadMmHA4Hl112GR9//DH3338/q1ev\nZv369Vx44YW0t7f3Wk9CQgIWy9G0Fhnp3EEnJycfs6zveydMmNDreWxsLPHx8QPezu/w4cPY7XaS\nkpKw2Wyexw033EB3dzeVlZWkpKTw9ttv09XVxdVXX01GRgYXX3wx+/fvH/mX5SW9zZ5SfRTXtTJ9\nwsivTO3JPTvkgZoWFuSm+GSd/pQeH8nCvBQKi+tZmJdCevzoj168NXPmTF588UW6urr48MMP+f73\nv8/FF1/M+++/z6ZNm3jjjTe44IILPO3b2nx7/UDfo4XW1laam5uZNGlSv+1TU1OJiIhgzZo1vXYo\nbu6dxaJFi3jzzTdpa2vjnXfe4fbbb+faa69l7dq1Po2/L03uSvVgdxjK6to4b5Zv7n2alRKDCJTW\nt/pkff4mIqy6ZRG1LZ2kx0f65OhluGw2G+ecc44nCbp7zlFRR68WLi4uZs2aNV7fVNsbb7/9Ns3N\nzZ6hmX/84x+ICAUFBf22P+ecc7Db7TQ0NHDeeecNuf6YmBguvfRStm7dyq9+9SufxT0QTe5K9VDV\n2E6n3THqGne3qAgrExOiKa0bP1epWixCRsLIp10Yic2bN3PnnXfyla98halTp1JfX8+vf/1rTjzx\nRBYtWkR2djZ33HEHP/vZz2hqauLee+8lK8u3U1zFxMRw8cUXc9ddd1FZWcldd93FFVdcwezZs/tt\nP2PGDFasWMHSpUu5++67KSgooL29nW3btrF7924ee+wxXnvtNZ544gkuv/xycnNzKS8v509/+lOv\ncwD+osldqR6Ka51li/lpoy+DdMtNjR03PfdAyczMZOLEifziF7+goqKC5ORkzj77bH79618TFRXF\nSy+9xK233spVV11FdnY2P/7xj3n//fd9OhXA0qVLSUhI4Oabb6a5uZnLLruMRx55ZND3PPzww0yf\nPp1HH32Ue+65h8TERGbPns3NN98MOE/eigg/+tGPqK6uJiMjg0suuYRf/vKXPot7IGICdK1RQUGB\n2bBhQ0C2rdRAnllXwg9e2sKHd59NziinH3C7/bnP+GRfLZ/88FyfrK+vHTt2MGvWLL+sO1zk5+dz\n1VVX8dvf/jbQoQCD/5+KSKExpv+xoh60WkapHorrWrFZhcnJw79v6kByUmKpamyno9vus3UqNRRN\n7kr1UFLbSnZKLFaL704k5qTGYgxUHGkfurFSPqJj7kr1UFzXQq6PhmPcclKcRwGlda2jnkZY+cfB\ngwcDHYLPac9dKRdjDMW1raOeMKwvd+VNSZ2eVFVjR5O7Ui5HWrtoau/2ec99YkI0kVaLXytmAlUY\noXzPV/+XmtyVcil29azzfFgGCc668ayUGMr8VOtus9l8frWmCpy2tjZsNtuo16PJXSkXd427r4dl\nALJTYvzWc58wYQLl5eW0trZqD34cM8bQ2tpKeXn5MfPcjISeUFXKpcQzj7vvk3tuaiyvb+l/AqrR\nSkxMBKCiooKurvEx+6Tqn81mY+LEiZ7/09HQ5K6US3FdKxMTo4i2WX2+7pzUWOpbu2hq7yIhevSH\n3H257/KjlJsOyyjlUlrX6pdeOzgvZHJuQ8fG1djQ5K6US1l9mycJ+1pOqqvWXeeYUWNEk7tSQLfd\nQVVjO9kpvpt2oKejPXdN7mpsaHJXCqhsaMfuMGT7qeeeHGsjISpCk7saM5rcleLocIm/eu4iQnZq\nLKX1OuauxoYmd6VwjrcDfuu5g3OOGe25q7EyZHIXkRwRWS0i20Vkm4h8t582IiIPisheEdksIgv8\nE65S/lFW34ZFYFJytN+2keO6aYdeaKTGgjc9927gDmPMbGARcKuI9L3v1IXA8a7HcmDw25coFWTK\n6luZlBSDzeq/g9mclBjauxzUNHf4bRtKuQ35m2yMqTTGbHT93ATsAPrevHAJ8FfjtBZIFpH+bxmu\nVBAqq28jy0/j7W7u2SG11l2NhWF1U0QkH5gPfNrnpSygtMfzMo7dASgVtMrr2/x2MtXNXQ5ZprXu\nagx4ndxFJB54EfieMaZxJBsTkeUiskFENtTU1IxkFUr5XJfdQWVDm19PpsLRk7VlWjGjxoBXyV1E\nbDgT+9PGmJf6aVIO5PR4nu1a1osxZqUxpsAYU5CRkTGSeJXyucoj7TiM/8og3WIiraTFRWrPXY0J\nb6plBHgc2GGM+Z8Bmr0CfM1VNbMIaDDG+GcKPKV8rMzPNe49ZafEaM9djQlvZoU8DVgGbBGRz1zL\nfgTkAhhj/gi8DlwE7AVagRt9H6pS/uFOtv6aV6an7JRYdlSOaFRTqWEZMrkbYz4CBr0VvHEW7t7q\nq6CUGktl9a1YBDKT/Ffj7paVEsPbOw7hcBgslkH/rJQaFb1CVYW9svo2v9e4u2WnxNDZ7eBwi9a6\nK//S5K7C3ljUuLu5x/V13F35myZ3FfbK6lvH5GQqQFaylkOqsaHJXYW1zm7nPO5jcTIV8BwhlGty\nV36myV2FtcqGtjGpcXeLj4ogJdamte7K7zS5q7A2FlP99pWlte5qDGhyV2FtLC9gcstOjqX8iCZ3\n5V+a3FVYK6tvw2oRJo1Bjbub8ypVnddd+ZcmdxXWyurbyEyMJmIMatzdsl3zute2dI7ZNlX40eSu\nwlppXeuY1bi7ZenskGoMaHJXYa2svm3MyiDdsrUcUo0BTe4qbHV02znU1E5O6lj33N1XqWo5pPIf\nTe4qbJXXt2HM2MwG2VNitI3E6AgdllF+pcldhS3PVL+pY5vcwVlXr+WQyp80uauwVRqAGnc3dzmk\nUv6iyV2FrdK6NmxWYWLi2NW4u7mvUtVad+UvmtxV2CqrbyUrOQZrAG6akZ0SS2unnSOtXWO+bRUe\nNLmrsFVa3zamc8r0pPO6K3/T5K7CVlld65iXQbplJWs5pPIvTe4qLLV2dlPb0hmwnnuOXqWq/EyT\nuwpLgSyDBEiMiSAhKkJ77spvNLmrsFRaF7gySAARITs1llLtuSs/0eSuwpKn5x6gYRnntmMoqdOe\nu/KPIZO7iDwhItUisnWA188SkQYR+cz1uMf3YSrlW6V1rUTbLKTHRwYshtzUWJ3XXfmNNz33vwAX\nDNHmQ2PMSa7H/aMPSyn/Kq1vJTslFpGxr3F3y0mNpb3LQU1zR8BiUKFryORujPkAqBuDWJQaM6V1\nbeQEaLzdzV2GWVqn4+7K93w15r5YRDaLyBsiMsdH61TKb8rqWwNWKeN2tBxSx92V70X4YB0bgVxj\nTLOIXAS8DBzfX0MRWQ4sB8jNzfXBppUavoa2LhrbuwN6MhXw1NiX6klV5Qej7rkbYxqNMc2un18H\nbCKSPkDblcaYAmNMQUZGxmg3rdSIBLoM0i0m0kp6fJQOyyi/GHVyF5FMcZ2VEpGTXeusHe16lfKX\nQF/A1FNOqpZDKv8YclhGRFYBZwHpIlIG3AvYAIwxfwSuAr4pIt1AG7DUaG2XCmLuMe5AD8u4Y9hY\nUh/oMFQIGjK5G2OuGeL1h4CHfBaRUn5WWtdKQlQEiTG+OOU0OjmpMby2pZJuu4MIq15TqHxHf5tU\n2CmrbyM7NbA17m65qbHYHYbKhvZAh6JCjCZ3FXZK61sDXuPulqMVM8pPNLmrsOJwGIprW8lPjwt0\nKMDRk7qlWuuufEyTuworVY3tdHQ7yEsL/MlUgElJ0VgtouWQyuc0uauwcrC2BYD8tODouUdYLUxK\nitZySOVzmtxVWCmudSbRYOm5g3PcXYdllK9pcldh5eDhFiIjLExOCo4TquAsh9RhGeVrmtxVWDlY\n20JuaiwWS+DLIN1yUmI53NxBW6c90KGoEKLJXYWV4tpW8oNoSAYgN01nh1S+p8ldhQ2Hw3CwtiVo\nTqa6eWaH1OSufEiTuwob1U0dtHc5yAuSGnc3vWmH8gdN7ipsHC2DDK5hmYz4KKJtFk8lj1K+oMld\nhY2Dh4Orxt1NRMhPi6PYtfNRyhc0uauwcbC2FZtVmJwcPGWQbvlpcRzQ5K58SJO7ChsHD7eQkxqL\nNYjKIN2mZMRRUttKt90R6FBUiNDkrsLGwdoWpgTZkIzblLQ4uh2G8iN6UlX5hiZ3FRaMcc4GmRek\nyd09S+X+wzo0o3xDk7sKC9VNHbR12clPD65KGbcpruR+UJO78hFN7iosuJNmsPbc0+MjiY+K0OSu\nfEaTuwoL7hryYB1zFxGmpMfpsIzyGU3uKiwcqG0hwiJMTo4OdCgDyk+P81xopdRoaXJXYaG41lkG\nGWEN3l/5KelxlNe30dGts0Oq0Qve33SlfOjA4eCbDbKvKemxOIzeLFv5xpDJXUSeEJFqEdk6wOsi\nIg+KyF4R2SwiC3wfplIj53AYDh5uCZqbYg9kSno84NwRKTVa3vTc/wJcMMjrFwLHux7LgUdGH5ZS\nvlNW30Zbl50ZExMCHcqg3Cd7DxxuDnAkKhQMmdyNMR8AdYM0WQL81TitBZJFZJKvAlRqtHYdagJg\nemZwJ/ekWBspsTbtuSufiPDBOrKA0h7Py1zLKn2wbjUObKto4K8fF5MSF8kNi/PJTAquipTdruR+\n/IT4AEcytCnpcdpzVz7hi+TuNRFZjnPohtzc3LHctPKDotIj/O6d3azeVUNspJX2LjuPf7Sfy0/K\n4jvnHk9OanCcwNxV1URWcgwJ0bZAhzKk/PQ4Pt5bG+gwVAjwRbVMOZDT43m2a9kxjDErjTEFxpiC\njIwMH2xaBcrW8ga+/MdPKCpr4M7zp/PJD8/lP3edzXWn5PHq5gqWrlzLkdbOQIcJOHvuM4J8SMZt\nanocVY3ttHZ2BzoUNc75Irm/AnzNVTWzCGgwxuiQTAhr7ujmtr9vJDUuknduP5PbzjmepBgbOamx\n/PSyOTz3jVOpbmrnjueKMMYENNYuu4N9Nc1MD/KTqW75njlmdNxdjY43pZCrgE+AGSJSJiI3i8gK\nEVnhavI6sB/YCzwKfMtv0aqAM8bwk39soaSulQevmU9qXOQxbeZlJ/Pji2bx7s5qHv1wfwCiPOrg\n4Ra67IYZmcE/3g49JhDTK1XVKA055m6MuWaI1w1wq88iUkHt+cIyXv6sgtvPm87JU1IHbHf94nzW\nHazj12/uYmFeCgvzBm7rT55KmfHSc/eUQ2pyV6OjV6gqr9U0dfDTV7ax+Lg0bj172qBtRYQHrpxH\nVnIMdz2/GbsjMMMzu6uasAgclzE+eu5xURFMTIxif40mdzU6mtyV1x77aD/tXXZ+fvlcr25Vlxht\n40cXzWT/4RZeLaoYgwiPtetQE/npcUTbrAHZ/khMn5jgKd9UaqQ0uSuv1Ld08tQnxVwybzJTh9EL\nPn92JjMzE3jwvT0B6b3vPtQc9Fem9jXDldwDdbSjQoMmd+WVP398kJZOO7edM/hwTF8Wi/Ddc49n\nf00L/9o8tr339i47B2tbxs14u9uMzAQ6uh16UlWNiiZ3NaTG9i7+vOYAF8zJHFGi/OIcZ+/99++O\nbe99b3UzxsDMcVLj7jZrUiLgvPhKqZHS5K6G9LdPimlq7x52r93NYhG+E4Deuzs5BvucMn1NmxCP\nRWBnZWOgQ1HjmCZ3Nai2TjuPfbifs2dkMDcracTruWBOJjMmJvDQe3vH7MKm3YeaiIywkBck0yB4\nK9pmJT89jp3ac1ejoMldDepfmyuob+3iG2ceN6r1WCzCLWdMZU91M2v3DzbJqO/sOtTEtIz4oL77\n0kBmZiZ4avSVGonx91uvxtQz60uZmhHHKYNcsOStS+ZNIjnWxt/WHhx9YF7YXTV+5pTpa2ZmIiV1\nrbR06BwzamQ0uasB7apqorC4nmtPzkVk6Lr2oUTbrFxdkMO/tx3iUGO7DyIc2OHmDioa2pk1aXwm\n9xmZCRiD1rurEdPkrga0al0JkVYLX1qQ7bN1XndKLg5j+PunJT5bZ38Ki+sBWJiX4tft+Iu7wkcr\nZtRIaXJX/WrvsvPSxjIumJvZ7+RgI5WXFseZ0zNYta6ELrvDZ+vtq7C4nkirhTmTR34SOJByUmKJ\njbTqSVXOE8FvAAAZXElEQVQ1YprcVb9e31JJY3s315zs+5uqfO3UPKqbOnh7+yGfr9utsLieE7KT\nxtW0Az1ZLML0iQnsrNJySDUymtxVv1atK2FKehyLpvp+Nsczp08gOyWGv35y0OfrBudRx5ayBgrG\n6ZCM28zMBHZVNQV8Tnw1PmlyV8fYW93E+oP1LP1cjk9OpPZltQjXnZLH2v117K32/f1Ct5Y30Gl3\nsGCcJ/cZmQnUt3ZR09QR6FDUOKTJXR3jhcJyrBbx6YnUvr5ckI3NKqxa5/sTqxvG+clUN3cZp467\nq5HQ5K56sTsM/9hUxlnTM8hIiPLbdtLjozh/TiYvbiyjvcvu03UXFtczJT2O9Hj/xT8WZmY655jR\ncXc1EprcVS9r9h7mUGMHVy70X6/d7bqTcznS2sUbW313y11jDBuL61mQO7577QCpcZFMSIjSnrsa\nEU3uqpcXN5aRFGPj3FkT/L6tU49LY0p6nE9r3g8cbqG2pZOC/PGf3AHmZiWxuawh0GGocUiTu/Jo\nau/i39uquPTESURF+L+EUES45uQc1h+s99mVmO6Ll8Z7pYzbwrwU9lY3c6S1M9ChqHFGk7vyeH1L\nJe1dDq7044nUvq5amEOk1eKz3nthcT2J0RHj5p6pQ3EPL20sqQ9wJGq80eSuPF4oLGNqRhwn5SSP\n2TZT4yK5YK7zxKovJsnaUFzPwrwULF7c43U8OCknGatFPEckSnlLk7sCoLi2hfUH67lyQbZfatsH\nc/3ifJrau3mhsGxU66luamdvdTMF+b6/8CpQYiKtzJmcyIaDmtzV8HiV3EXkAhHZJSJ7ReQH/bx+\nlog0iMhnrsc9vg9V+dMLhWVYBL60IGvMt70wL4UFuck8sebAqG7D99pmZ9XN+bMn+iq0oLAgN4Wi\nsiN+nYtHhZ4hk7uIWIGHgQuB2cA1IjK7n6YfGmNOcj3u93Gcyo/sDsMLhWWcMT2DSUkxAYnh66dP\npbi2dVTzzbxSVMGsSYkcP85uiD2UgvwU2rsc7NDb7qlh8KbnfjKw1xiz3xjTCTwDLPFvWGosfbT3\nMJUN7Xx5YU7AYjh/9kSyU2J4/KP9I3p/SW0rm0qOcNmJk30cWeC5r7TVcXc1HN4k9yygtMfzMtey\nvhaLyGYReUNE5vgkOjUmnt9QSnKsjS/M9n9t+0AirBZuOm0K6w/W81npkWG//1XXjbcvPXGSr0ML\nuElJMUxOivZMq6CUN3x1QnUjkGuMmQf8AXi5v0YislxENojIhpqaGh9tWo3GkdZO3tp2iMtPyhqT\n2vbBXP25HBKiI3jsw+H33l/5rIKCvBSyU8bXzbC9tTA/lY2a3NUweJPcy4Gex+vZrmUexphGY0yz\n6+fXAZuIpPddkTFmpTGmwBhTkJGRMYqwla/887MKOu0Ori4I3JCMW3xUBNeenMsbW6vYM4yLmnZW\nNbLrUBOXnRR6QzJuC3OTqWxop+JIW6BDUeOEN8l9PXC8iEwRkUhgKfBKzwYikimu+jkROdm13lpf\nB6t877kNpczNSmT25MRAhwLAN848jvioCH766jav5zF/5bMKrBbhohNCb0jGbWGes7xTx92Vt4ZM\n7saYbuA24N/ADuA5Y8w2EVkhIitcza4CtopIEfAgsNToHQaC3tbyBrZVNAZFr90tNS6SO784gzV7\na3lja9WQ7Y0xvLq5gtOmpY/7WSAHM2tSAjE2qyZ35bUIbxq5hlpe77Psjz1+fgh4yLehKX97am0x\n0TYLS04c+9r2wVx7ci6rPi3h5//azlkzMoiNHPjX9IXCMkrr2rj9vOljGOHYi7BamJ+bzNr9ekCs\nvKNXqIapI62dvPxZOVfMzyIp1hbocHqxWoT7l8yhoqGd/1u9b8B2pXWt3Pfqdk6Zkhp0Oyh/OHvG\nBHZWNVFa1xroUNQ4oMk9TD27vpT2LgfXL84PdCj9KshP5Uvzs1j5wX5W76w+5nWHw3Dn80UA/PbL\nJ4bMXDKDOc915e1bfryxuAodmtzDkN1h+NvaYk6Zkuq5208w+skls5meGc/NT67nb2uLe732xJoD\nfHqgjnsunU1OamiWP/aVnx7HjIkJvL196HMRSnk15q5Cy7s7DlFW38aPL5oV6FAGlRoXybPLT+Xb\nqzbx3y9vZXtFI4nREWyraOTTA7WcN3siXx6DO0YFk/NmT+T/3t9LfUsnKXGRgQ5HBTHtuYehJz85\nyOSkaM9hfjCLi4pg5bKFfO3UPFatK+HPaw5ypK2Tqwty+PWV88Z8BstAO3/ORBwG3u1nqEqpnrTn\nHmb2HGpizd5a7vriDCKs42PfHmG1cP+Sudx29jRS4iKxBUPcFZuG/57J80e92ROykshMjOatbVVc\nFWZHLWp4guCvRI2lR97fR7TNwtLPBU9tu7cmJEYHR2IPIBHhvNkT+WBPDW2d9kCHo4KY9tzDyN7q\nZl7+rJyvnz6VtBC+4CdoDbe3P0BP//w5E/nb2mI+2nt4XAytqcAI725QmHnw3T1E26x844ypgQ5F\njcIpU9JIiI7grW1aNaMGpsk9TOw+1MSrmyu4YXG+9trHucgIC1+YNZE3t1XR7IP7zqrQpMk9TPzu\nnd3ERUZwy+naaw8FXzs1j6b2bp5bXzp0YxWWdMw9DGyvaOT1LVV855xp/qmNDlDlSDibn5tCQV4K\nT6w5wNdOzRs3lU9q7OhvRIhzOAz3vrKVpBgbN39ee+2h5OunT6Wsvo1/b9PpCNSxNLmHuKc/LWb9\nwXp+cvGsoJsgTI3OebMnkpcWy6Mf7vd67nsVPnRYJoRVHGnjgTd2cvrx6XrBy1BGMrQUYFaLcPPn\np3DPP7dRWFxPQX5qoENSQUR77iHKGMNPXt6Kw8Avrzgh7C7TDxdXLcwmKcbGyg+Gf99ZFdo0uYeo\nFzeW897Oau784oywmTUxHMVGRnDjafm8tf0Qq3fpfDPqKE3uIWjdgTp+9NIWTp6Syg1BOl+78p0V\nZx7H9InxfP+FzRxp7Qx0OCpI6Jh7iNlX08wtf91AdkoMf/rqQqxhcBOLfo3DMfSRirZZ+Z+rT+Ly\nh9fw3//cxh+u0TJTpT33kFLT1MENf15HhEX4y40n63zfYWRuVhLf+8LxvFpUwStFFYEORwUB7bmH\niL3VTSz/ayE1TR08s/xUctNGMc4eRr1ef3A4oKG9i+QYG2N5HnvFmcfxzo5qfvyPLWQlR7MwT6tn\nwpkm9xDw721V3PFcEdE2C3+7+RROykn26n0Oh6G2pZP0+MhBq2n8kqx8NENifxwOQ0Pr2CdX57bh\nRy9vYUdlI7MmJfLLy0/AYhmbhB9htfDQtfNZ9vg6rnvsUx66ZgFf0Fkjw5YOy4xj9S2d3P/qdr7x\nt0KOy4jj1W9/ns95WevscBiueXQtp/7qXZauXIvD0f9FMO5kdcOf1/HDf2zB4fAuNocD6lu7cF9b\n0/e5v7g/13DjHXyd3sfe0N7FjspG7A7DjspGGtq7+v0O/fV9ZKfE8sKKU5k+MYFvPFXIM+tK9AKn\nMOVVz11ELgB+D1iBx4wxD/R5XVyvXwS0AjcYYzb6OFbl0tLRzZ/XHOBP/9lPc2c3156Syz2XzCba\nZu3/Df30khtau2gvKWSmMbSXCA37raT0cwVrr2RV0UBxXQv5aXGe3mfPHqkxzp8ToyP48ctbPb3X\nny+Zy0/+ufWY3uywDNDT79sjdn+unsm1v8812Dr6vubuic/MTOTuC2aQGhuJSP/vS46xMWtSouez\nJkZHUFzXwo6KBuwGdlQ2UtfSwa/e2MGeQ83MmpzIzy+bS+mRVpJibKTGRR2NYYRHN2nxUay6ZREr\nnirkBy9t4R+byvn+hTNZkJsyvPWpcW3I5C4iVuBh4DygDFgvIq8YY7b3aHYhcLzrcQrwiOtf5SNt\nnXb+s7ua17ZU8e6OQ7R22jlv9kTuPH8GMzIThr2+vkkoOab/BOhpV9FAdGQE33v2M0+Cht6JDww7\nq5o4fkI8e6qbPQm2tL61V2+2vq0Ti4hXQxTeJl53TN58rr47pL7rgKOv99y5bato4MYn1jE7K4mf\nL5nLj1/ews7KRmZOSuQXl59AU0c38ZFWlp8+haQYG8mxkZ6dXGSElY5uOzMzE7j/1a3sr20DYFtF\nI9c+/iltXc67Ks2ZlMAvrphHU0e35zP3/Q76+04cDqht6vAMscVFRfDEDZ/j75+W8If39vCl//uY\nL8yayOXzJ3PG9AwSo3UqilDnTc/9ZGCvMWY/gIg8AywBeib3JcBfjfP4b62IJIvIJGNMpc8jDlHG\nGFo77TS2d3GosYPKI21UNLSzs7KRLeUNnmSZGhfJkpOyuLogm/nD6In1TQgi8MvLTxhyHNjdrriu\nhe8+swmHgR0VDTS0dwF4Et/OKue/ADurmoiJtNLWaSfaZiUnJcazg5iWEcdv3tzJzqqmfpOpu/c/\nUOLtOX5tjDlmCCQl1jbo5+q7Q7jrizN673haO/l/b+1y7bASuOv8GcycGM+2yibn+4Ht5Q0U17aw\nraIRcCboO577jIN1rRhjcBjnjb0fvma+Z91tDmfy3lXZSHefURJ3YgfYVtnEHS8UcfBwCzMnJfLz\nS+fww5e3sOdQE7MmJ/V7FATOz/RCRSUL81JYdcsiLBbBZrVw/eJ8rlqYzeMfHeDPaw7wzo5DRFiE\ngvwUZk9KYtqEeKZmxDEhIYq0uCgSYyL0auYQ4U1yzwJ6ThpdxrG98v7aZAE+T+4f7zvM797e4+vV\njprh6F+sMWBwJmzjfu76o7c7DHaHocvhoLPbQUe3g44uOy2ddk9y7CktLpK5WUl8YdZEFk1NY9HU\n1GFP7zrQST6LhSGHLMDZLjc1lpjICFo6uomOjCAxOgKLiKeXHGm1eBKYRfDc37O9s5umDjs/XzKX\n77+0mT2HmnB/zB0VDdS1dPDbt3d7kikIO6tciff8Gcck76Ro29GjhUmJzMxM8OwokmNsQ5647Dsm\nLtCrpy9Cj556Izc9uZ7IiN7fd5TNyu/f2dVr2b7DLb2et3R009TexbQJ8eyqavIs75vYoyxgsVo9\nCT7KKuyrbgZgW3kDd764mX01zZ7vq+9RUM+dbLcjjcLiempbOslIOHpDlrioCL5z7vHcevY0NpXU\n886OatbsPcyqdSW9dizgnK8mNtJKtM1KjM2KzSpEWCxEWAWrRRARLAKC836uzn+Pvt+5RA3lspMm\n89VFeX7dxphWy4jIcmA5QG5u7sjWgQTphTnS+5dcnLE6e8lH/yCsFgsRFsFqFaIiLERFWImKsBAf\nFUFiTAQJ0TYmJESRmRTN5KQYkmNto+5J1SbN5oWKSrodaWyrEO5Imt3rj7+/qpm+y+qaOtjQWYXd\nGKydQl3SHNLiIvmvZbM43NzOxX9Yc3SDBmIihPZuw0nZSSRNLaC2pYt/Hqqi25HRq90t79jZXJ6G\n3aSxtcL5xdldcd6eeSLRuQ4Ki+tZmJtC8nGf43BzZ6/Psub752CxCOnxkdjthqtXfkJR6REK8lM9\nPdieko0hOtfOhuJ65k1OojPjBL63bB4WETISonA4DJbJXWwubXDurg3Q96LPTuDw4N95XKSFvBNO\n4zdzT+OqP33CppIj/Te0g9UIsyYn8MAVc7ns4Y/puY/fXgMOMwGABTnJ5M09lehPrL2+E4DoXDsR\nxfUszEshPb7/6xusFqEgP9UzwZjDYahoaOPA4RYON3dQ19JFXUsHrZ122rvstHXa6XIYuu0Ouu0G\nh6uD4jDG1YExvU4I63lb741FDvMmuZcDOT2eZ7uWDbcNxpiVwEqAgoKCEf0qnHpcGqcelzaSt4at\n9PhIFualOBNCnz9+d3WJ+7VVtywCOGZZenwkBT3WkRpr45pH17LhYB2xkcf+GrW5uqibShtYuvJT\nVt1yiieGE7IS2Vzu7H1uLm/gxJxkNpc1sCA3GUTY6NpGRkIUT998Cntrmpk+MR4ROeazTEiMQkRw\nOAxfXnk0iW44WHdMDxacO9qnbz6FL//pYzaVHmHRA6sBODk/hb9/fRHXPraWTaUNAMRHRdDa0U1c\ndARN7d1YcA7LeOOd/zoDi+us8YsrFrO94giXPPSx5/XZmXHsrGrBgfNobmdlExMSYyjIT2XdgToA\nCnKTsFidiXxedhIvrDgVi8XCqlsWHbMz7m/ZUCwWITslluwUnXsoFHmT3NcDx4vIFJwJeylwbZ82\nrwC3ucbjTwEadLw9OLh74H//+inUtXYd88df29JJYXE93Q7jOaQHjlmWFhfJg0vnIwIZCVEcbna+\nz26gaYj7eBaW1LOnutkTQ1qcjWse/dSToHvGZgyeJGUMXPf4p712MhaL9JvIals62VzW4NnmiTnJ\n/fZgHQ7D3ppminq0BdhwsJ491U0UFtd7lrV12Xn9u6d7ThBf9OCHuEffLALzspOwWSwUltTTd0RN\nxNLr6GdOVjIFuckUlhxhYV4yzy4/ldqWTm5btanXzuyZWxZR09Th+Z57fh/uz2qxyDE7rf6WqfA2\nZHI3xnSLyG3Av3GWQj5hjNkmIitcr/8ReB1nGeRenKWQN/ovZOWt/nrlfXt1A/XqF/bTS+/bk1+Y\nl+LsuUc5e7YA8ZFW2rodxNgsNHc4x3NjI61c/OCHvYZK+iZod2JyJzWAw80dx+xkMhKi+k1k7iOL\nDT16uX0/q/v72HCwjliblebOo+PNDuDeV7axIDeZ9cVHPN/BjMwERIQZmQl8Lj+VwuJ6FuSm8NC1\n8z3J93BzB8v/uoHPeu4w5Nijn+dWLO71mSckRvNMn+9BBCYmRR9dTY/vQ6nhkEBd4FBQUGA2bNgQ\nkG2Hi5qmDk791bt0OwwRFuGTH57bb6IYasz9cHNnv+txt0mNtXG4uRMR5wngutYu57KWTuqaO7jk\nDx9hNwwaQ3+MMSxdeTRBPrP82J3TUJ9joO/DIkdPfPe09gfnIBbnacGMhKhe6xls/Xa7g6v++DFF\npQ0U5Kfw0LULWPzAe0N+90oNl4gUGmMKhmqn0w+EsMHG2nsa6jB/oPX0bNOzt+lZlhjNhIQoClw9\n3sFi6I9I/0MwAxlqaKLn51iQlwLGsOFgfa9xdIvF2aMe7vqtVgsvfvM0T6yAV9+9Uv6iPfcQ5+38\nMf5cj69i8IWesbiHVG79+0bPuPez3zh2OMcX2wr051ahw9ueuyZ3FfY0CavxRIdllPKSVpqoUKSz\nQiqlVAjS5K6UUiFIk7tSSoUgTe5KKRWCNLkrpVQI0uSulFIhSJO7UkqFoIBdxCQiNUBxQDbev3SG\nnKU7oII9Pgj+GIM9PtAYfSHY44PRxZhnjMkYqlHAknuwEZEN3lz1FSjBHh8Ef4zBHh9ojL4Q7PHB\n2MSowzJKKRWCNLkrpVQI0uR+1MpABzCEYI8Pgj/GYI8PNEZfCPb4YAxi1DF3pZQKQdpzV0qpEKTJ\nvQ8RuUNEjIikBzqWvkTk/4nIThHZLCL/EJHkQMcEICIXiMguEdkrIj8IdDx9iUiOiKwWke0isk1E\nvhvomPojIlYR2SQi/wp0LP0RkWQRecH1O7hDRE4NdEx9ich/uf6Pt4rIKhHp/7ZaYxvTEyJSLSJb\neyxLFZG3RWSP698UX29Xk3sPIpIDnA+UBDqWAbwNzDXGzAN2Az8McDyIiBV4GLgQmA1cIyKzAxvV\nMbqBO4wxs4FFwK1BGCPAd4EdgQ5iEL8H3jTGzAROJMhiFZEs4DtAgTFmLmAFlgY2KgD+AlzQZ9kP\ngHeNMccD77qe+5Qm997+F7ibY++bHBSMMW8ZY7pdT9cC2YGMx+VkYK8xZr8xphN4BlgS4Jh6McZU\nGmM2un5uwpmUsgIbVW8ikg1cDDwW6Fj6IyJJwBnA4wDGmE5jzJHARtWvCCBGRCKAWKAiwPFgjPkA\nqOuzeAnwpOvnJ4HLfb1dTe4uIrIEKDfGFAU6Fi/dBLwR6CBwJsnSHs/LCLLE2ZOI5APzgU8DG8kx\nfoezY+EYqmGATAFqgD+7ho4eE5G4QAfVkzGmHPgtziPvSqDBGPNWYKMa0ERjTKXr5ypgoq83EFbJ\nXUTecY3F9X0sAX4E3BPkMbrb/BjnUMPTgYt0/BGReOBF4HvGmMZAx+MmIpcA1caYwkDHMogIYAHw\niDFmPtCCH4YSRsM1br0E545oMhAnIl8NbFRDM86SRZ+PFoTVPVSNMV/ob7mInIDzF6LIdYPkbGCj\niJxsjKkawxAHjNFNRG4ALgHONcFRx1oO5PR4nu1aFlRExIYzsT9tjHkp0PH0cRpwmYhcBEQDiSLy\nlDEmmBJTGVBmjHEf8bxAkCV34AvAAWNMDYCIvAQsBp4KaFT9OyQik4wxlSIyCaj29QbCquc+EGPM\nFmPMBGNMvjEmH+cv8oKxTuxDEZELcB66X2aMaQ10PC7rgeNFZIqIROI8gfVKgGPqRZx77MeBHcaY\n/wl0PH0ZY35ojMl2/e4tBd4LssSO62+hVERmuBadC2wPYEj9KQEWiUis6//8XILspG8PrwDXu36+\nHvinrzcQVj33EPAQEAW87TrCWGuMWRHIgIwx3SJyG/BvnNUJTxhjtgUypn6cBiwDtojIZ65lPzLG\nvB7AmMajbwNPu3bi+4EbAxxPL8aYT0XkBWAjzmHLTQTB1aoisgo4C0gXkTLgXuAB4DkRuRnn7LhX\n+3y7wXFkr5RSypd0WEYppUKQJnellApBmtyVUioEaXJXSqkQpMldKaVCkCZ3pZQKQZrclVIqBGly\nV0qpEPT/AZUnG5w5PnMJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb7b8bf8c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def func(x):\n",
    "    return np.exp(-x ** 2) + 3 * np.exp(-(x - 3) ** 2)\n",
    "x = np.linspace(-5, 10, 100)\n",
    "rv = Gaussian(mu=np.array([2.]), var=np.array([2.]))\n",
    "\n",
    "samples = metropolis(func, Gaussian(mu=np.zeros(1), var=np.ones(1)), n=100, downsample=10)\n",
    "plt.plot(x, func(x), label=r\"$\\tilde{p}(z)$\")\n",
    "plt.hist(samples, normed=True, alpha=0.2)\n",
    "plt.scatter(samples, np.random.normal(scale=.03, size=(100, 1)), s=5, label=\"samples\")\n",
    "plt.legend(fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11.2.2Metrppolis-Hasting算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 设定状态转移矩阵T,设定转移次数阈值$n_1$,需要的样本数量$n_2$\n",
    "    \n",
    "* 从任意简单概率分布采样得到初始状态矩阵$p(z^{(0)})$\n",
    "    \n",
    "* $m=0$到$m=n_1+n_2-1$\n",
    "\n",
    "     * 从条件概率分布$q_k(z^{(m+1)}|z^{(m)})$采样得到样本$z^{(m+1)}$\n",
    "        \n",
    "     * 从均匀分布采样$u-uniform(0,1)$\n",
    "     \n",
    "     *如果$$u<min(1,\\frac{\\tilde p(z^*)q_k(z^{(\\tau)}|z^*)}{\\tilde p(z^{(\\tau)})q_k(z^*|z^{(\\tau)})})$$则接受转移$z^{(\\tau)}\\rightarrow z^*$\n",
    "     * 否则不接受转移，$\\tau=max(\\tau-1,0)$\n",
    "* 样本集$(z_{n_1},...z_{n_1+n_2+1})$即为平稳分布对应的样本集"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(x):\n",
    "    return np.exp(-x ** 2) + 3 * np.exp(-(x - 3) ** 2)\n",
    "x = np.linspace(-5, 10, 100)\n",
    "#提议分布\n",
    "rv = Gaussian(mu=np.array([2.]), var=np.array([2.]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9//HXZyb7RhKSkJCVRcIioJAiQrUWqxU33Nqi\nlrq1lmq369baPqq39qrt/d17+6utV4tLa6sVRdT6K2jdl6oIBAir7GQjISF7Mllnvr8/ZiYkIcsk\nnMlMks/z8ZgHmTNnzvnMJLznO9/zPd8jxhiUUkqNLrZAF6CUUsp6Gu5KKTUKabgrpdQopOGulFKj\nkIa7UkqNQhruSik1Cmm4K6XUKKThrpRSo5CGu1JKjUIhgdpxUlKSycnJCdTulVJqRMrPzz9ujEke\naL2AhXtOTg6bN28O1O6VUmpEEpFCX9bTbhmllBqFNNyVUmoU0nBXSqlRSMNdKaVGoQHDXUQiRGSj\niBSIyC4R+WUv64iIPCIiB0Rku4jM80+5SimlfOHLaJlWYIkxplFEQoF/icjrxpgNXdZZCpzmuZ0F\nPOb5VymlVAAMGO7GfammRs/dUM+t5+WblgF/8ay7QUTiRSTNGFNmabVKqV7V19dTUVFBe3t7oEtR\npyA0NJSUlBTi4uJOeVs+jXMXETuQD0wFHjXGfNZjlXSguMv9Es8yDXc1apXVNfPxgSquOjMdm00C\nVkd9fT3Hjh0jPT2dyMhIRAJXixo6YwzNzc2UlpYCnHLA+xTuxhgncIaIxAOviMjpxpidg92ZiNwK\n3AqQlZU12KcrFTQKimv59l82U9nQSlSYnYtnpwWsloqKCtLT04mKigpYDerUiQhRUVGkp6dz9OjR\nUw73QY2WMcbUAu8BF/V4qBTI7HI/w7Os5/NXGWPyjDF5yckDnj2rVFBav6OMb6z6lPAQGznjo/jt\nW/twugJ3ofn29nYiIyMDtn9lrcjISEu613wZLZPsabEjIpHABcDnPVZ7DfiWZ9TMQqBO+9vVaPTm\nrnJue24LM9PiePX2xdz11Vz2VzSybkdg/9y1K2b0sOp36UvLPQ14T0S2A5uAt4wx/xCRlSKy0rPO\neuAQcAB4ArjNkuqUCjJrt5SQNi6Cv31nIUkx4Vx8ehq5E2L5v28HtvWuVE++jJbZDpzZy/LHu/xs\ngNutLU2p4NLW4eLjA1VcNnciEaF2AGw24d8uOI2Vz27h79tKuWpeRoCrVMpNz1BVykdbimpobO3g\nS9O6Hy+6cGYqM9Pi+N07++lwugJUnVLdabgr5aMP9lUSYhMWTR3fbbnNJtz25SkUVjnYdKQmQNWN\nHjfccAMJCQn87ne/4/nnn+fGG2/s9rgxhjPOOINnnnnG521+//vf55ZbbrG40uCm4a6Ujz7YW8m8\n7ATiIkJPemzxlCQAthZruJ+Kmpoa/vnPf/Lkk0/y7LPPcvvtt3PzzTd3W+fFF1+kurqa6667zuft\n3nXXXTz33HMcOHDA6pKDloa7Uj6oqG9hd1k95+X2PoQ3ITqMyUnRbCmsHebKRpeEhATKy8u5+uqr\n2bRpE9XV1Zx77rnd1nnkkUdYsWIFoaEnf8j2JScnhy9+8Ys89thjVpcctDTclfLBh/uPA5zU397V\nmVkJbC2qwT2+QA1FYmIiv/jFL7j++usZP348iYmJfPe736W1tRWAAwcO8Mknn3DNNdd0e94LL7yA\niPR6q6+vB+Dqq6/mueeew+UaG8dFNNyV8sH7eytIjg1nZlrfZw2emRVPVVMbxdXNw1jZ6FFUVERN\nTQ2//e1vSUxMZM2aNfz4xz9m1apVPProowC88847REdHM3fu3G7PPffcc/n00087b+vWrSMxMZGL\nL76480zPRYsWcezYMXbs2DHsry0QAnYNVaVGCqfL8NH+43xlxoR+TzCZl5UAuEfVZI3XqQAGq6Cg\nAHAf/Pz1r38NwJIlS/jwww954403uOOOO8jPz2fGjBnYbN3bpWlpaaSluaeAqK+v5/zzz+f0009n\nzZo1nevMmjULu93Oxo0bT/pwGI003JUaQEFJLXXN7Xypj/52r9zUWKLC7GwtquGKM9OHqbre/fL/\n7WL30fqA7HvmxDjuv2zWoJ9XUFBAVFQUd999d7fl06ZNY+PGjQCUl5eTlJTU5zaampq4+OKLsdls\n/OMf/+g2305ISAjx8fGUl5cPuraRSMNdqQF8tO84InDO1L5DBcBuE+ZmxLOlSA+qDkVBQQHnnnsu\n48d3H2paXl7OxIkTAWhpaelzgrSWlhYuv/xympqaePfdd4mNjT1pnfDwcFpaWqwvPghpuCs1gF1H\n65iUFE1CdNiA687LjuePHxyiuc1JZJh9GKrr3VBazoFWUFDAOeec021ZQ0MD7777Lr/8pfsCcImJ\nib22vNvb27n66qspKyvjgw8+ICEhodd91NbWkpiYaH3xQUgPqCo1gH3HGpieenIrsDfzshLocBl2\nlNb5uarRpampiYMHD3Lw4MFuyx955BHAfWITQG5uLocPH+62jtPp5Nprr2Xv3r28/fbb9DXjbGVl\nJQ6Hg2nTpvnhFQQfDXel+tHc5qSw2sG0Cb6F+xmZ8YD7oKry3Y4dO3C5XBw8eJAHHniA9957j3vv\nvZdf/vKXrFq1ivh49/u6ePFiioqKqKys7Hzu9773PdavX89DDz1EUVERGzZsYMOGDSeNitm8eTMi\nwqJFi4b1tQWKdsso1Y/9FQ0YA7k+hvv4mHByxkexVcN9UAoKCoiLi2PdunXcdNNNPPjgg+Tm5rJm\nzRqWLVvWud55551HYmIib7zxBitWrMAYw+rVq2lubuYb3/hGt21edtllvPbaa53333jjDb70pS+d\n1Kc/Wmm4K9WPveUNgHskjK/mZSXw0YHjGGN0nnUfFRQUMHv2bObMmUN+fn6f64WFhfHNb36T1atX\ns2LFim4nKfXH6XSydu3aziGWY4F2yyjVj73lDYSF2MgeH+3zc87MiqeyoZWSGj2ZyVcFBQU+jz2/\n++67ee+999i3b5/P21+zZg2RkZEsX758qCWOOBruSvVj77EGTkuJwT6IC2DPyXD3D+8K0DjzkcYY\nw44dO5gzZ45P62dkZPD0009TVub71a+MMTz11FOEhIydzoqx80qVGoJ9xxpYPMD49p4mJ7tb+YeO\nN/qjpFHH166VrgbbAr/22msHtf5ooC13pfpQ62jjWH2rzwdTvWIjQpkQF86hyiY/VabUwDTcleqD\n92DqtEEcTPWanBTDwUptuavA0XBXqg/7jrnD3dcTmLqakhLNocomnf5XBYyGu1J9+Ly8gdiIEFLj\nIgb93MlJMdQ1t1PV1OaHypQamIa7Un3wTjswlLHqnQdVtd9dBYiGu1K9MMawt7zB52kHepqSHAPA\nIe13VwEyYLiLSKaIvCciu0Vkl4j8qJd1zhOROhHZ5rnd559ylRoe5fUt1Ld0DOrM1K7S4yMJD7Hp\nQVUVML6Mc+8A7jTGbBGRWCBfRN4yxuzusd5HxphLrS9RqeHXOe3AEFvuNpswKSlau2VUwAzYcjfG\nlBljtnh+bgD2AIG9zIxSftY5DHKI4Q7urhltuatAGVSfu4jkAGcCn/Xy8CIR2S4ir4vIyLtSgFJd\n7K9oJDk23KcLdPRlcnI0xTXNtHW4LKxMKd/4HO4iEgOsBX5sjOl5rvAWIMsYMwf4PfBqH9u4VUQ2\ni8jmrvMxKxVsCquamDSIycJ6MyU5BqfLUFStXTOjwfvvv4+IsHPnzkCX4hOfwl1EQnEH+3PGmJd7\nPm6MqTfGNHp+Xg+EishJE3IYY1YZY/KMMXl9XS1FqWBQWOUge3zv1+r0lXc45IEKDXc1/HwZLSPA\nU8AeY8z/9LFOqmc9RGSBZ7tVVhaq1HBxtHVQ0dBqQbh7hkPqBGIqAHxpuS8GVgBLugx1vFhEVorI\nSs861wA7RaQAeARYbvS8azVCFVU7AMg6xW6ZmPAQJsSFc1Bb7j7ZtWsXF110EYmJiURHRzNjxgwe\nffRRANatW8cFF1xASkoKcXFxLFy4kDfffLPb8//93/+dpKQkPvvsM/Ly8oiMjOSLX/wihw8fpqKi\ngiuuuIKYmBhmzJjBu+++2+25OTk53HXXXfzqV78iNTWVmJgYrr/+eurq+r8Wrsvl4te//jVTp04l\nPDycadOm8cwzz3Rb51//+hfnnHMOcXFxxMXFccYZZ7BmzRoL3rH+DTgU0hjzL6DfU/SMMX8A/mBV\nUUoFUmGVO9yzE0+t5Q7uaQi05e6byy67jBkzZvDss88SHh7O3r17O6cCPnz4MBdffDF33nkndrud\n119/naVLl/Lhhx+yePHizm04HA5uvfVW7rnnHqKjo/nhD3/IihUrCA8PZ+nSpdx2223853/+J1/7\n2tcoLi4mKurE7/j5559n6tSpPPHEE5SVlXHPPffw7W9/u98g/sEPfsAzzzzDfffdx7x583jrrbe4\n+eabGT9+PJdeein19fVceumlLFu2jPvuu69z7vra2lr/vZFexpiA3ObPn2+UCkarPjhosn/yD1PT\n1HrK2/r5K9vN7PvfMC6Xy4LKerd7925Lt+d0ukxFfYtfa+6psrLSAGb79u0Drut0Ok17e7u58MIL\nzU033dS5/P777zeAef/99zuXPfroowYwv/zlLzuX7dq1ywBm/fr1ncuys7NNQkKCaWho6Fz27LPP\nGhHpfH/fe+89A5gdO3YYY4zZv3+/ERHz5z//uVt9K1asMHl5ecYYYzZt2mQAU19fP5i3o9/fKbDZ\n+JCxOv2AUj0UVjcxLjKU+KihD4P0mpwUQ31Lx4iZQMzlMlz7xAbOfvgdlq/agMs1PL2riYmJZGZm\nsnLlSl544QUqKiq6PV5SUsINN9xAeno6ISEhhIaG8uabb550qb2wsDDOOeeczvtTp04FYMmSJSct\nKy0t7fbcCy64gJiYmM77V155JcYYNm3a1GvN77zzDjabjSuvvJKOjo7O2/nnn8+2bdtwOp1MmTKF\nmJgYrrvuOv7+978PT4vdQ8NdqR6sGCnjNdImEKtqaiO/sIYOlyG/sGbYPpRsNhtvvvkmqamp3Hzz\nzaSmpnLOOeewdetWXC4Xl19+OZ988gkPPPAA7733Hps2bWLp0qW0tLR0205sbCw224lYCwtzf0DH\nx8eftKznc1NSUrrdj4qKIiYmps/L+R0/fhyn08m4ceMIDQ3tvN144410dHRQVlZGQkICb731Fu3t\n7Xz9618nOTmZSy65hEOHDg39zfKRXmZPqR6Kqh3MTh9nybYmJ7lbgkeON7FgUqIl2/SnpJgw5mcn\nkF9Yw/zsBJJiTv3bi6+mT5/O2rVraW9v56OPPuInP/kJl1xyCe+//z5bt27l9ddf56KLLupcv7nZ\n2guQ9/y24HA4aGxsJC0trdf1ExMTCQkJ4eOPP+72geLl/bBYuHAhb7zxBs3Nzbz99tvccccdXHfd\ndWzYsMHS+nvScFeqi3ani5KaZi6d0/t/6MFKi4/AbhOKaxyWbM/fRITnv7OQqqY2kmLChjTd8akK\nDQ1lyZIlnSHobTmHh4d3rlNYWMjHH3/s80W1ffHWW2/R2NjY2TXzyiuvICLk5eX1uv6SJUtwOp3U\n1dVxwQUXDLj9yMhILrvsMnbu3MnDDz9sWd190XBXqoujtc04XYbsxFMbBukVareRNi6ic3jlSGCz\nCcmx4QOvaKHt27dz11138Y1vfIPJkydTU1PDb37zG+bOncvChQvJyMjgzjvv5Fe/+hUNDQ3cf//9\npKdbO8VVZGQkl1xyCXfffTdlZWXcfffdXHnllcycObPX9XNzc1m5ciXLly/nnnvuIS8vj5aWFnbt\n2sW+fft48sknWbduHU8//TRXXHEFWVlZlJaW8sc//rHbMQB/0XBXqovOYZAW9bkDZCZEUTyCwj0Q\nUlNTmTBhAg8++CBHjx4lPj6eL3/5y/zmN78hPDycl19+mdtvv51rrrmGjIwMfv7zn/P+++9bOhXA\n8uXLiY2N5ZZbbqGxsZHLL7+cxx57rN/nPProo0ybNo0nnniC++67j7i4OGbOnMktt9wCuA/eigg/\n+9nPqKioIDk5mUsvvZSHHnrIsrr7IiZA5xrl5eWZzZs3B2TfSvXlrxsK+cWrO9lw7/mkjhv85fV6\n85OXtvPu3go2/fwrlmyvpz179jBjxgy/bHusyMnJ4ZprruG//uu/Al0K0P/vVETyjTG99xV1oaNl\nlOqiqKqJ8BAbKRZ2S2QmRlLZ0Epzm9OybSo1EA13pboorHKQlRiFzWbdgcRMz5muJSPkoKoaHbTP\nXakurBzj7uUN9+IaB6edwsU/lP8cOXIk0CVYTlvuSnkYYyiqdpB9ihOG9ZSZ4A73oiptuavho+Gu\nlEdlQyvN7U7LW+5JMWFEhtoprrH2pJuuAjUwQlnPqt+lhrtSHoXeqX4tmA2yKxEhMzHSb8MhQ0ND\nLT9bUwVOc3MzoaGhp7wdDXelPE6Mcbe2WwbcXTP+OpEpJSWF0tJSHA6HtuBHMGMMDoeD0tLSk+a5\nGQo9oKqUR1FVE3abkB4fafm2MxOj+OxwNcYYy0/pj4uLA+Do0aO0t7dbum01vEJDQ5kwYULn7/RU\naLgr5VFY7WBifARhIdZ/oc1MjKKxtYNaRzsJ0dZPxuW9yo9SXtoto5RHcbXD8v52r8wE97eBkTTH\njBrZNNyV8iipaSYj3k/h3mWsu1LDQcNdKaCl3UlFQysZCdb3t0OXcK/WUS1qeGi4KwWU1rpDNyPR\nP+EeEx5CYnSYdsuoYaPhrhTuLhmAjAT/dMuAu99d55dRw0XDXSlOTOqV6c9wT9R53dXwGTDcRSRT\nRN4Tkd0isktEftTLOiIij4jIARHZLiLz/FOuUv5RXN1MqF0sneq3p8zEKEo9V3pSyt98abl3AHca\nY2YCC4HbRaTndaeWAqd5brcC/V++RKkgU1LjID0+0tKpfnvKTIii3Wkor2/x2z6U8how3I0xZcaY\nLZ6fG4A9QM+LFy4D/mLcNgDxImLNFYaVGgYlNc1+7W+HE3PW6OyQajgMqs9dRHKAM4HPejyUDhR3\nuV/CyR8ASgUtd7j7Z6SMV6ZnJI6OdVfDwedwF5EYYC3wY2NM/VB2JiK3ishmEdlcWVk5lE0oZbmW\ndifHG1s7x6L7S9q4SEROjMxRyp98CncRCcUd7M8ZY17uZZVSILPL/QzPsm6MMauMMXnGmLzk5OSh\n1KuU5bwjZfzdcg8LsZEaF0GphrsaBr6MlhHgKWCPMeZ/+ljtNeBbnlEzC4E6Y0yZhXUq5TfFnWPc\n/Rvu3n3oWHc1HHyZFXIxsALYISLbPMt+BmQBGGMeB9YDFwMHAAdwk/WlKuUfw3ECk1dGQhQbD1f7\nfT9KDRjuxph/Af2ODzPuKwTcblVRSg2nkhoHYSE2kmP8N8bdKyMhktcKWuhwugix6zmEyn/0r0uN\neSXVzWT4eYy7V3p8JE6XoaxOx7or/9JwV2NeSY2D9GHob4cTXT/eicqU8hcNdzXmDccJTF7eg7Y6\nHFL5m4a7GtMcbR1UNbUNy0gZgLT4CM9Ydx0xo/xLw12NaaXDOAwSIDzETkpsuLbcld9puKsxzTsV\ngL/PTu0qIyFKT2RSfqfhrsa0kmFuuXv3VVKr3TLKvzTc1ZhWUtNM+DCNcffKSIikrNY91l0pf9Fw\nV2Oadxike5aN4ZGREEWHy3CsoXXY9qnGHg13NaYVVw/fMEiv9HjPcEi95J7yIw13NaYV1zjIHMb+\ndjjRv68nMil/0nBXY1ZDSzu1jvZhHSkDMDFeT2RS/qfhrsYsb7hmDnO3TESod6y7dsso/9FwV2NW\ncfXwXKSjN+553bXlrvxHw12NWd6LdAx3twxAekKU9rkrv9JwV2NWcbWD6DA7CVGhw77vjIRIjtY2\n43SZYd+3Ghs03NWYVVLjIDMxaljHuHtlJETS7jRUNOi87so/NNzVmDWcU/325N2v9rsrf9FwV2OS\nMYbiakdADqZClxOZdMSM8hMNdzUm1TjaaWpzBuRgKnS5aEe1ttyVf2i4qzHJOwxyuM9O9YoItZMc\nG9455bBSVtNwV2NSIOZx7ykzIZJibbkrP9FwV2NSSQDHuHtlJkZpy135zYDhLiJPi0iFiOzs4/Hz\nRKRORLZ5bvdZX6ZS1iqudpAQFUpMeEjAashKjOJobTPtOq+78gNfWu5/Bi4aYJ2PjDFneG4PnHpZ\nSvlXcQCHQXplJkThMlBWq2PdlfUGDHdjzIdA9TDUotSwKal2kJkYmIOpXhme/WvXjPIHq/rcF4nI\ndhF5XURmWbRNpfzC5TKU1DQP+2yQPXn3X6wX7VB+YEWH4xYgyxjTKCIXA68Cp/W2oojcCtwKkJWV\nZcGulRq8ysZW2pwuMgJ4MBUgbVwEdptQpOGu/OCUW+7GmHpjTKPn5/VAqIgk9bHuKmNMnjEmLzk5\n+VR3rdSQBHqMu1eI3UZ6fGTn7JRKWemUw11EUsUz85KILPBss+pUt6uUvwTDGHevzMRI7ZZRfjFg\nt4yIPA+cBySJSAlwPxAKYIx5HLgG+J6IdADNwHJjjM5jqoKW98Qh7/wugZSZEMXbe44Fugw1Cg0Y\n7saYawd4/A/AHyyrSCk/K652kBIbTkSoPdClkJkYxfHGNppaO4gO4Jh7NfroGapqzCmpaQ6KLhk4\n0TWkU/8qq2m4qzGnqNoR8IOpXt46tN9dWU3DXY0pLe1OjtY1k5MUHehSgBMtdz2RSVlNw12NKSU1\nDoyBnPHBEe7jo8OIDLXrWHdlOQ13NaYcPu4O0WBpuYsIWYlROvWvspyGuxpTCquaAMgZHxwHVME9\n1l0vt6espuGuxpTDx5uIjwolPios0KV0ykiIorjagZ4eoqyk4a7GlMIqB9lB0t/ulZkYRVObk+qm\ntkCXokYRDXc1phw+3hRUXTLQZTikjnVXFtJwV2NGa4dnGGSQtdyzxuvUv8p6Gu5qzHD3a0NOUrC1\n3HWsu7KehrsaM454h0EGWcs9OjyE8dFhFB7XcFfW0XBXY8aRzmGQwRXuAJOSojnsqU8pK2i4qzHj\nSFUT4yJDSYgOnmGQXjlJ0Rw5ruGurKPhrsaMwipH0I2U8ZqUFE1FQyuNrR2BLkWNEhruasw4fLwp\n6Ma4e03yTIegrXdlFQ13NSa0djg5Whs8s0H25A33wxruyiIa7mpMKK5uxmWCa06ZrrwHebXlrqyi\n4a7GhM4Jw4K05R4ZZidtXIS23JVlNNzVmOANzWAcBuk1KSmaQxruyiIa7mpMKKxyEBcRQkJUaKBL\n6VNOUnTnWHylTpWGuxoTjlQ1kZMUjYgEupQ+TU6KptbRTo3ODqksoOGuxoRgHgbp5e0y0jNVlRUG\nDHcReVpEKkRkZx+Pi4g8IiIHRGS7iMyzvkylhs7R1kFpbTNTk2MCXUq/JiV7wr1Sw12dOl9a7n8G\nLurn8aXAaZ7brcBjp16WUtY5UNGIMZCbGtzhnpkQhU3QfndliQHD3RjzIVDdzyrLgL8Ytw1AvIik\nWVWgGhnqW9ppaXcGuoxefV7eAMC0CbEBrqR/YSE2MhOjdMSMskSIBdtIB4q73C/xLCuzYNsqyB05\n3sQfPzzI2vxSosPtfOvsHL51djbjY8IDXVqnfeUNhIfYgr7PHdzDIfVEJmUFK8LdZyJyK+6uG7Ky\nsoZz18piLpfh56/u4IVNxYTYbVyTl0FFfSu/e2c/j39wkLu/msu3z5kc6DIB2HusgdMmxGC3Be9I\nGa+c8dFsPFyNMSaoR/ao4GdFuJcCmV3uZ3iWncQYswpYBZCXl6eXeh/BHvvgIM9vLOaGs7O5fclU\nUmIjADhQ0cBD6z/nP9btYWpKDOflpgS4Uth3rIHFU5MCXYZPJidH42hzUtnQSkpcRKDLUSOYFUMh\nXwO+5Rk1sxCoM8Zol8wotvFwNf/95l4unzuRf798VmewA0xNieXR6+YxPTWWf3thG2V1gb3oc62j\njWP1reQGeX+7l3c4pPa7q1Ply1DI54FPgVwRKRGRW0RkpYis9KyyHjgEHACeAG7zW7Uq4Gqa2vjR\n6q1kJUbx4JWn99p1EBlm59Hr59HW4eIHf9tKu9MVgErd9h1rBGBa6sgId536V1llwG4ZY8y1Azxu\ngNstq0gFLWMMd79UQFVjGy/ftojYiL5P5Z+SHMNDV83mR6u38du39nHPRdOHsdIT9h5zj5QZKS33\nifGRhNltOoGYOmV6hqry2ScHq3h7TwV3fzWX09PHDbj+sjPSuWpeOk9+dDhg3TP7yhuIjQghbdzI\n6L+224RJSdEcqGgMdClqhNNwVz575J39TIgL51uLsn1+zh0XTMNlDI+9f9CPlfVt77EGcifEjqiR\nJ7mpsZ1j85UaKg135ZONh6v57HA13z13CuEhdp+fl5EQxdfyMlm9sXjYW+/GGPaWN4yY/nav3NRY\nSmubqW9pD3QpagTTcFc++f27+0mKCePaBYM/P+H2L08JSOu9oqGVuub2EdPf7jUjzV3vPm29q1Og\n4a4GtK24lo/2H+c750wmMsz3VrtXoFrve0fItAM95abGAbBHw12dAg13NaDfv7Of+KhQrl/oe197\nT4Fove/zjpQZYd0yE8dFEBsRwt7y+kCXokYwDXfVrwMVDbzzeQU3L55ETPjQT2jOSIji6nkZvLi5\nmDrH8PQl7y1vIDk2nMTosGHZn1VEhOmpsXxepi13NXQa7qpff/usmFC7cN1Zpz4X0A2Lcmhpd7Em\nv3jglS2wzzNSZiTKTY1l77EG3KeRKDV4Gu6qTy3tTtZuKeGrs1JJsmCWx5kT48jLTuCvGwpxufwb\nWh1OF/uONY64/nav6alxNLR0cLSuJdClqBFKw1316Y2d5dQ1t3PdEEbI9GXF2dkUVjn46MBxy7bZ\nm8/LG2hudzI3c+CTrYLRdM9xgs/LtN9dDY2Gu+rT3zYWkTM+ioWTx1u2zaWnp5EUE8ZfPz1i2TZ7\nk19YA8D87AS/7sdfvGPz9WQmNVQa7qpXByoa2Hi4muULsrBZOA96WIiN5V/I4p3PKyiudli23Z42\nF9aQGhdBenyk3/bhT3ERoaTHR2q4qyHTcFe9en6j+0DqNfMzLN/2dWdlIcBznxVZvm2v/CPVzM9J\nGFHTDvQy4bxVAAAYnElEQVQ0PTVWh0OqIdNwVyfxHki90KIDqT1NjI/kgpkTeHFzMa0d1l939Wht\nM0frWsgboV0yXtPTYjlY2eSX90iNfhru6iRv7T5GraOda7/gv0shfnNhNtVNbbyxs9zybXv72/Oy\nEy3f9nDKTY3D6TIcrNDpf9Xgabirk6zdUsLEcREsmmLdgdSeFk9JIisxyi9dM/mFNUSG2jvnaBmp\nvCNm9h7Trhk1eBruqptj9S18uK+Sq+ZlWHogtSebTbh2QRYbD1dzoMLag4abC6s5IzOeEPvI/vOe\nlBRNmN2mZ6qqIRnZf/3Kcq9uLcVl4Kp56X7f19fyMgi1C3/7zLozVptaO9hT1kBezsjubwcItduY\nmhLDbh3rroZAw111MsbwUn4J87LimZwc4/f9JcWE89VZqbyUX0xLuzUHDbcV1+J0mRE7vr2nM7Pi\n2Vrkfk1KDYaGu+q0o7SO/RWNXDM/c9j2ef1Z2dS3dLBue5kl29t8pAYRmDdKwn1+dgKNrR2dM1wq\n5SsNd9VpbX4JYSE2LpmTNmz7XDg5kcnJ0Tz3WaEl29tcWE3uhFji+rl490jiHfGz2TMCSClfabgr\nAFo7nPy94ChfnZXKuMjhC0YR4fqzstlSVMv2ktpT2pbTZdhWVDtqumQAMhMjSYoJZ4uGuxokDXcF\nwLt7Kqh1tHP1MBxI7enreRnEhofwxEeHT2k7nx2qoqG1g0VTkiyqLPBEhLzshM6x+0r5yqdwF5GL\nRGSviBwQkZ/28vh5IlInIts8t/usL1X504ubi0mNi+Cc05KHfd+xEaEsX5DJ+h1llNYO/TJ8rxUc\nJTrMzvkzUiysLvDmZydQVO2gokGn/1W+GzDcRcQOPAosBWYC14rIzF5W/cgYc4bn9oDFdSo/Kq9r\n4YN9lVwzPwO7H8e29+fGxZMA+NO/htZ6b+1wsn5HGV+dlUpE6OCv8xrM5nuGdWrXjBoMX1ruC4AD\nxphDxpg2YDWwzL9lqeG0dksJLuMedx4o6fGRXDw7jdWbimloGfxl+D7YW0l9SweXnzHRD9UF1qyJ\ncYSF2Nh8RMNd+c6XcE8Hup5lUuJZ1tMiEdkuIq+LyKzeNiQit4rIZhHZXFlZOYRyldWMMby4uZiF\nkxPJHh8d0Fq+c84kGls7eGHT4E9qeq3gKInRYSyeOnr6273CQ+zMzRhHfpGGu/KdVQdUtwBZxpg5\nwO+BV3tbyRizyhiTZ4zJS04e/r5ddbKNh6sprHLw9bzhG9velzkZ8SyYlMifPj5CW4fL5+c1tnbw\n9p5jXDI7jdARPuVAX+ZlJ7CztM6yk73U6OfL/4RSoOv//AzPsk7GmHpjTKPn5/VAqIiMvibUKPTC\n5mJiwkNYevrwjW3vz/e/PJXS2mae/tj3vve3dpfT0u5i2SjskvHKy06k3WnYUVoX6FLUCOFLuG8C\nThORSSISBiwHXuu6goikiueqCCKywLPdKquLVdZqaGln/Y4yLps7kciw4DgIee60ZL4yYwKPvLOf\nch8vDv3atqOkx0cyL2v0jG/vaV5WPID2uyufDRjuxpgO4PvAP4E9wIvGmF0islJEVnpWuwbYKSIF\nwCPAcmOMToYR5F7ddpSWdhff+ELgu2S6uv+ymXS4DA+t3zPguiU1Dj7af5zL5k706yyWgTY+JpzJ\nSdHkF1YHuhQ1QoT4spKnq2V9j2WPd/n5D8AfrC1N+ZMxhr98coTZ6eOYmzEu0OV0k5kYxfe+NIXf\nvbOfaxdkcXYf88q7XIa712wnPMTG9Wf578IiweLsKeN5ZWspLe3OUTfcU1lvdB59UgP65GAV+ysa\nuWFRTlBeZ/R7500hIyGS+1/bSXNb7wcRn/74MJ8equK+y2aSmRg1zBUOvwtnpeJoc/LJweOBLkWN\nABruY9SfPj7C+OgwLh3GScIGIyLUzgPLZrG/opHlT2ygsqG12+P7jjXwn//cy1dmpATFSJ/hcPbk\n8cSGh/DmrmOBLkWNABruY1BxtYN3Pj/GtQuygvrr/ZLpE3j8m/PZW17Plf/7MQcqGqhsaOX9vRX8\naPU2YsNDePiqOUH5zcMfwkJsfCk3mbf3HNP53dWAfOpzV6PLXzcUYhPh+oXB30/91VmpvHDr2dzy\nzGYu/O2HeDMt1C48dv18kmPDA1vgMLtwVir/2F7G1qIa8nJG9gXAlX9puI8xjrYOVm8s4qJZqaSN\niwx0OT6ZmxnPK7ct4q8bCkmNi2DmxDhmpMUN69TEAzq6dfDPmXjmoJ9yXm4yoXbhzd3HNNxVvzTc\nx5i1+SXUt3Rww6KcQJcyKJmJUfzs4hn+39FQQnoYxUWEcvaUJN7cVc69S6ePmS4pNXja5z6GtLQ7\nefS9g+RlJ/CFUXAB6bHqwpkTOFLl4EBFY6BLUUFMw30MeX5jEeX1Ldxx4TRt8Y1gF8ycAMCbu3XU\njOqbhvsY0dzm5H/fP8jCyYmj6kpFY9GEuAjOyIxn/Y4y9ERw1RcN9zHi2Q2FVDa0cscFuYEuRVng\n6vkZ7DparxfOVn3ScB8Dmlo7ePyDg5xzWhILJukIi9HgmnkZJESFsurDQ4EuRQUpDfcxYNWHh6hq\nauPfLpgW6FKURSLD7KxYmM3be45xqFIPrKqTabiPcvuONfC/7x/g8rkTR/WUuGPRirNzCLXbeGqI\n151Vo5uG+yjmdBnueWk7sRGh3H9Zb9c0VyNZcmw4V52Zzkv5JVQ1tg78BDWm6ElMo9ifPj7MtuJa\nfrf8DMbH9HGa/lBP2hnC2ZXKet8+ZxKrNxXz1w2F/Pgr2u2mTtBwH6WKqhz815t7WTI9hcvnjt7L\nz/UpyM80BSz5YJ2aEsuS6Sn8+ZMjXH9W9piba0f1TbtlRqHG1g5WPptPiM3Gf1xxup6wNMrdu3Q6\njjYn9768Xce9q07ach9l2p0ubntuC3uPNfDUDXlMjPfT5GDDNFGWGthpE2K556u5/Me6PazJLxkz\n89ur/mnLfRQxxvCLV3fy4b5KHrzidM7LTQl0SWqY3Lx4EmdNSuSB/7eb4mpHoMtRQUDDfZRwugy/\nfv1zVm8q5gdLprJ8QfDP1a6sY7MJ//31uQDc8eI2Wtp7vzShGjs03EeBuuZ2vv3MJv744SG+uTCL\nO3w8WcnlMtQ42umtm9blotfH+lo+lgTre5CREMWDV57O5sIavvnkZ9Q62gJdkgog7XMf4XYdreP2\n57ZQWtvMf1xxOtefleXTAVSXy3DtExtoKcpnRlocD10xG5vN+xjc+8oOPi+vZ3pqLPdcNJ3EqDCM\ngZ+9uoM9ZfUnPWckc4d1GyKQEBVGf2+fy2Xde+ByQV1LO/GRoZ377G3ZYCw7Ix27TbjjhQK+9vin\nPHPzAv8dd1FBTcN9hCqudvDbt/fxytZSxkeH8/x3Fg7qyjxVTW3kF9Yw3Rj2lNVT19JOQpT7ykY1\njjZ2Ha0DYNfRem58eiO5qbHcu3Q6e8rqcbpOfs6phlKgeD/IvK931sQ4Hr5yDtD766lrae/zPRjs\nfnt+SLiM4Scvb+dAReMpfXBcOmciidFh3PqXfK549GN+ctF0rjjTHfpq7PDpT0dELhKRvSJyQER+\n2svjIiKPeB7fLiLzrC9VtTtdfLCvkrvXFLDkv99n3fYyvnPOZN76t3OZl5VAZUOrz0PhkmLCmJ+d\ngN0mzEiLI77LJet6hrMBPi9v4KH1e5ieGovdJkzv8hyXC+59uYAbnvqMn64twOWy5vW6XGZQr2ng\n7Z3cnVLX0s7n5fWd9/ccrafa0crPXt3BjX/ayI9f3EZlQ0vnc+IjQ5mRFofdJkxNiWFcREiv+3C6\nDFWNbVQ3tXU+t+v+e35IHKlq5CdrC9hb3tDtg6O/2vuzaEoSL373bFLiwrlzTQEX/+4j/rmrnHan\nRb8cFfQGbLmLiB14FLgAKAE2ichrxpjdXVZbCpzmuZ0FPOb5Vw2RMYb65g72lNezs7SOHaV1fLCv\nklpHOzHhIXwtL5MfLJlK2rjIzi6W/MIa5mcn8Px3FmIboJUmIjz/nYXUHbJ3tk67tr5nTYxjT1k9\nYSE2WtrdgXCgsomc8VEYlwGXi+ONrdhtQrvTxa6yBgB2lTVQ5WglubczYvsZPtmz5d+1ZRuRNd+n\n1+Tr9rq2iuMjQ5meFseuUnfL3QU8vH4PByqbcLoMBysauenPm5mZGsvdS6djF+FXl8/ipy9vZ195\nA3e/VMCvr5pDfWsHGMN//nMvn5fVExEaQlNbB+D+NvDAZbO455UdHKpoJDc1ll9fNZvpqbF8XtZA\nRIidH63eRtfcnpoSQ1xECDWOduIiQvj5qzu71I5PLfqZwGtXRbNun53/s6mV7/41n9jwEM6dlsx5\nucnMnBjHlOQYIkLtPr+vauTwpVtmAXDAGHMIQERWA8uAruG+DPiLcTexNohIvIikGWPKrC64uNrB\npwerrN7sSQyDay12bVEZz32D8fzrDmuXy+A04HS5aHca2jpctHa4aO1w0tjSQUNLB/Ut7Ryrb6Gs\nrgVH24kRDxPiwjlvWjIXz07j3GnJ3f5DnuhiOUhLkVB3yO5TV4ENunWrdA2/B6+YTUNrB7Hhdn7y\n8g4OVDQSFmLjYGUT4A7xm/68CYDwkO6h29jcflK499dt01vwdm3ZthTl+/yafNnenrJ6ahxt2GxC\nfGQoDy47nTtf3MrB4+4hhHuPNZIzPpojVU2d29xd3sBNf3K/3twJMew91ti57rVPbKClo3uL2Bvs\n4O7auu6pz2j1rPN5eQM/WbuDouomnMZ0W9e7/d9cNbsz0E9LiWHfsQZcBnYfraOwuomc8dG+dX8Z\n4eNdByitGc9pKTHMy0rg3b0VrNvh/q8pAhPHRZISF8746DASosKIDg8hItROZKid0BAhxCaE2GzY\nbYJN3A0DERC8/54wlC45Yex1F+WmxjI3M96v+/Al3NOB4i73Szi5Vd7bOulAt3AXkVuBWwGysoY2\nVG97SR33rN0+pOcGozC7jbAQGzHhIcRFhhAbEcq0CbGcOy2ZieMimZoSw6z0OFJiI/rchreLJb9w\nCvOzEoif8oVB/y+ramjlpaNldLjGs7MUbgyZSm5GLCLCb24/k01Hqrn2yc96f/KJ3gOiw+xkz17c\nrWnZ85vFc7ecRU1zO0kxYYhIt33vOircOW4mSelhRGQ53c8Z5GuqrG9hTelRnGY8O48Kd8TNYHxa\nGPb0dnYW1zIvI4Fvv9XO9pI68nISeWT5mbxWWdbt4/zhb5zN1x7/lJZeRhTuLAfocg5B+8nr9Pce\nAewuBxfJJ61mF+Gpm75Mjc3GmtJynGY828sgKjyVplYnduD8v9XxhUl2n77NeN9bpzEcPt7E376z\nkIevms3+ikb2VzRwoKKRw8ebqGpso7S2hZ2l9TjaOmhpd9GmXTh+s/JLU4Ii3C1jjFkFrALIy8sb\nUkfqkukpfPzTJZbW1ZfBtie6Zk/XVo3IiRaP3eZuCdltQpjdNqiuhr736+5iqWpq6wzMwfJ+QGw+\nUk1UeAiX/P5f5HmC+PqnPmNzYQ3R4XaaWvsePz07PY5Xb1uErUew7zvWwObCGpwuw+bCGr6+6lO2\nl9R1diGd+HByh7/3NQzlNblchh88vxWn56/L6TLc/lw+IjYKSuqYnTGODpdha7G7G2bzkWoMhlkT\nY9l5tKFzO8v+91PmpI9jW0ndgPuMCrXhaB9cEEaG2WnyfDOLChVyU+PYWlyH07jrf+6Ws4gKD6Gh\npYOocHtn15j3deUX1lDV1DbgXDJ9vbe5qbHkpsb2+9wOp4sOl3HfnC5cBlyeb6Bdv516DeXQSJCN\nJh02MeH+j15f9lAKdD2fOcOzbLDrWCIyzE56mA7t6slmk1OaNMobpvuONXDJ7/+F02XIL6zhQGUj\n+Z5gbm5zMjs9jt1lDcxJH8fW4trO57/+wy8yPS2uWwh7W+zeDwxHm5M5GeMoKK7FaboHVG9BPpTX\n5O2i6iq/qBYRweky7Cip6/YpPCdjHD9cvY3dZQ1EhdlobnNhAJeBHUfrOSMjjm0l9d22F2GnW4t+\nzcqzuez3H9NXvNttQkSIrTPMgc5uHLsI7931ZRBh0cPv4DSwpaiWg8ebOrvlmtuczM2MZ3tJHVFh\ndhytHZ1BPZBT+eAPsdsI0e74EcuX0TKbgNNEZJKIhAHLgdd6rPMa8C3PqJmFQJ0/+tuVf9ls7hZd\nXnYCITZhfnYC0ybEMN9z/ws5ibx622I23Hs+a793NgtyErALLMhJOCnY4UTQOg04WjtY94Mvsnbl\n2eTlJHZu3xtQ3iA/1UnOkmLCmJMxrtuyORnjOl/D/OyEzlFCZ2bF8/g357OlsAaXgdZ2F+t/+EUW\nTHLXl5edwNrvLWbdDxZ3217XYD8zM54ZaXF9fsX+QnY8n/50CQX3XcAZGeOwCeRlx3fWk5eTQEpc\nBCmx4d3el2kTYjp/D1/ISeSl757NhnvPZ9svLmDDz77C6lsX+vxeWfXeqpFFfBlmJiIXA/8XsANP\nG2MeFJGVAMaYx8X9V/MH4CLAAdxkjNnc3zbz8vLM5s39rqKGkctlOlt3xtCtpdf1sZ4t8/5ahMYY\nlq860dfuDaSBnneqnE4XX/vjp2wrquWMrHheWnk2IL2+PuCkGnu+fmMM3/jjp+QX1jAvOwHB/a1j\nbqZ72zabjfK6ZhY+/G5nDVGhNt676zxS4iI6X/PyVRvIL3Lv5289jjv09n76+31SI5OI5Btj8gZc\nL1BThGq4B4+hDKUczLYDEVCD2a8v6/b34Qd0fgBsOlLDrImx/P32xdjtJ/o0KhtaOfvhd+hwGUJs\nwqf3nq9zr6sh8TXc9QxV1dl90uHpZ/flQJ2vTvVYwHDs15d1u64jwknriwirbz27zw+J3g5sKuVP\nGu5Kg8ci/X1IWDGiSanB0HBXGjzDJFDfYtTYpOGuAA0epUabUTBhq1JKqZ403JVSahTScFdKqVFI\nw10ppUYhDXellBqFNNyVUmoU0nBXSqlRKGBzy4hIJVAYkJ33Lgk4Hugi+hHs9UHw1xjs9YHWaIVg\nrw9OrcZsY8zJV3rpIWDhHmxEZLMvk/EESrDXB8FfY7DXB1qjFYK9PhieGrVbRimlRiENd6WUGoU0\n3E9YFegCBhDs9UHw1xjs9YHWaIVgrw+GoUbtc1dKqVFIW+5KKTUKabj3ICJ3iogRkaRA19KTiPwf\nEflcRLaLyCsi0vtVmYeZiFwkIntF5ICI/DTQ9fQkIpki8p6I7BaRXSLyo0DX1BsRsYvIVhH5R6Br\n6Y2IxIvIS56/wT0icnaga+pJRP7N8zveKSLPi0hEENT0tIhUiMjOLssSReQtEdnv+TfB6v1quHch\nIpnAhUBRoGvpw1vA6caYOcA+4N4A14OI2IFHgaXATOBaEZkZ2KpO0gHcaYyZCSwEbg/CGgF+BOwJ\ndBH9+B3whjFmOjCXIKtVRNKBHwJ5xpjTATuwPLBVAfBn4KIey34KvGOMOQ14x3PfUhru3f0WuAcI\nygMRxpg3jTEdnrsbgIxA1uOxADhgjDlkjGkDVgPLAlxTN8aYMmPMFs/PDbhDKT2wVXUnIhnAJcCT\nga6lNyIyDjgXeArAGNNmjKkNbFW9CgEiRSQEiAKOBrgejDEfAtU9Fi8DnvH8/AxwhdX71XD3EJFl\nQKkxpiDQtfjoZuD1QBeBOySLu9wvIciCsysRyQHOBD4LbCUn+b+4GxauQBfSh0lAJfAnT9fRkyIS\nHeiiujLGlAL/hfubdxlQZ4x5M7BV9WmCMabM83M5MMHqHYypcBeRtz19cT1vy4CfAfcFeY3edX6O\nu6vhucBVOvKISAywFvixMaY+0PV4icilQIUxJj/QtfQjBJgHPGaMORNowg9dCafC02+9DPcH0UQg\nWkS+GdiqBmbcQxYt7y0YU9dQNcZ8pbflIjIb9x9Egefi0BnAFhFZYIwpH8YS+6zRS0RuBC4FzjfB\nMY61FMjscj/DsyyoiEgo7mB/zhjzcqDr6WExcLmIXAxEAHEi8qwxJpiCqQQoMcZ4v/G8RJCFO/AV\n4LAxphJARF4GFgHPBrSq3h0TkTRjTJmIpAEVVu9gTLXc+2KM2WGMSTHG5BhjcnD/Ic8b7mAfiIhc\nhPur++XGGEeg6/HYBJwmIpNEJAz3AazXAlxTN+L+xH4K2GOM+Z9A19OTMeZeY0yG529vOfBukAU7\nnv8LxSKS61l0PrA7gCX1pghYKCJRnt/5+QTZQd8uXgNu8Px8A/B3q3cwplruo8AfgHDgLc83jA3G\nmJWBLMgY0yEi3wf+iXt0wtPGmF2BrKkXi4EVwA4R2eZZ9jNjzPoA1jQS/QB4zvMhfgi4KcD1dGOM\n+UxEXgK24O623EoQnK0qIs8D5wFJIlIC3A/8GnhRRG7BPTvu1y3fb3B8s1dKKWUl7ZZRSqlRSMNd\nKaVGIQ13pZQahTTclVJqFNJwV0qpUUjDXSmlRiENd6WUGoU03JVSahT6/+ohH2xU0KZZAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb7b54ebd30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = metropolis_hastings(func, Gaussian(mu=np.ones(1), var=np.ones(1)), n=100, downsample=10)\n",
    "plt.plot(x, func(x), label=r\"$\\tilde{p}(z)$\")\n",
    "plt.hist(samples, normed=True, alpha=0.2)\n",
    "plt.scatter(samples, np.random.normal(scale=.03, size=(100, 1)), s=5, label=\"samples\")\n",
    "plt.legend(fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "提议分布的具体的选择会对算法的表现产⽣重要的影响。对于连续状态空间来说，⼀个常见\n",
    "的选择是⼀个以当前状态为中⼼的⾼斯分布，这会在确定分布的⽅差参数时需要进⾏⼀个重要\n",
    "的折中。如果⽅差过⼩，那么接受的转移的⽐例会很⾼，但是遍历状态空间的形式是⼀个缓慢\n",
    "的随机游⾛过程，导致较长的时间开销。然⽽，如果⽅差过⼤，那么拒绝率会很⾼，因为在我\n",
    "们考虑的这种复杂问题中，许多的步骤会到达$p(z)$很低的状态。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "rv=multivariate_gaussian.MultivariateGaussian(np.array([5,-1]),np.array([[1,0.5],[0.5,2]]))\n",
    "def p_ygivenx(x,m1,m2,s1,s2):\n",
    "    return(random.normalvariate(m2+rho*s2/s1*(x-m1),math.sqrt(1-rho**2)*s2))\n",
    "def p_xgiveny(y,m1,m2,s1,s2):\n",
    "    return(random.normalvariate(m1+rho*s1/s2*(y-m2),math.sqrt(1-rho**2)*s1))\n",
    "N=5000\n",
    "K=20\n",
    "x_res=[]\n",
    "y_res=[]\n",
    "z_res=[]\n",
    "m1=5\n",
    "m2=-1\n",
    "s1=1\n",
    "s2=2\n",
    "\n",
    "rho=0.5\n",
    "y=m2\n",
    "for i in range(N):\n",
    "    for j in range(K):\n",
    "        x=p_xgiveny(y,m1,m2,s1,s2)\n",
    "        y=p_ygivenx(x,m1,m2,s1,s2)\n",
    "        z=rv._pdf([x,y])\n",
    "        x_res.append(x)\n",
    "        y_res.append(y)\n",
    "        z_res.append(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  2.68210561e-05,   0.00000000e+00,   2.68210561e-05,\n",
       "          1.07284225e-04,   2.95031618e-04,   6.43705348e-04,\n",
       "          8.85094853e-04,   1.55562126e-03,   2.60164245e-03,\n",
       "          5.25692700e-03,   8.52909585e-03,   1.24181490e-02,\n",
       "          2.07058553e-02,   2.93422354e-02,   4.25381950e-02,\n",
       "          5.55732283e-02,   7.41602202e-02,   9.69044759e-02,\n",
       "          1.18495426e-01,   1.42178419e-01,   1.62347853e-01,\n",
       "          1.79727897e-01,   1.91287772e-01,   1.97215226e-01,\n",
       "          1.98851310e-01,   1.89678509e-01,   1.76589834e-01,\n",
       "          1.60604484e-01,   1.40381408e-01,   1.18066289e-01,\n",
       "          9.56707073e-02,   7.45893571e-02,   5.92745341e-02,\n",
       "          4.14921739e-02,   2.95836249e-02,   2.04912869e-02,\n",
       "          1.30618543e-02,   8.63638008e-03,   5.25692700e-03,\n",
       "          3.11124251e-03,   1.93111604e-03,   1.04602119e-03,\n",
       "          3.75494786e-04,   3.48673730e-04,   1.34105281e-04,\n",
       "          5.36421123e-05,   2.68210561e-05,   0.00000000e+00,\n",
       "          0.00000000e+00,   2.68210561e-05]),\n",
       " array([-9.94280777, -9.56996638, -9.19712498, -8.82428358, -8.45144219,\n",
       "        -8.07860079, -7.7057594 , -7.332918  , -6.96007661, -6.58723521,\n",
       "        -6.21439382, -5.84155242, -5.46871103, -5.09586963, -4.72302823,\n",
       "        -4.35018684, -3.97734544, -3.60450405, -3.23166265, -2.85882126,\n",
       "        -2.48597986, -2.11313847, -1.74029707, -1.36745568, -0.99461428,\n",
       "        -0.62177289, -0.24893149,  0.12390991,  0.4967513 ,  0.8695927 ,\n",
       "         1.24243409,  1.61527549,  1.98811688,  2.36095828,  2.73379967,\n",
       "         3.10664107,  3.47948246,  3.85232386,  4.22516526,  4.59800665,\n",
       "         4.97084805,  5.34368944,  5.71653084,  6.08937223,  6.46221363,\n",
       "         6.83505502,  7.20789642,  7.58073781,  7.95357921,  8.3264206 ,\n",
       "         8.699262  ]),\n",
       " <a list of 50 Patch objects>)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFNNJREFUeJzt3W9sXfd93/H3Z3SdB5q7rJHiePpTK6jQzAPizOPkoNBS\nG409yRiqpO1WpUXTNgkEAVG7PAg2FQECFXnkDRuGAkoELRPaDk21AY06IWWixsFUP0jdiQ4c23Ks\nhFFcSIITyU6QLEgRW8h3D+6Rd8xcmofkvbwkz/sFEDx/fr/LL88lz4fnd/4wVYUkqX/+3qQLkCRN\nhgEgST1lAEhSTxkAktRTBoAk9ZQBIEk9ZQBIUk8ZAJLUUwaAJPXULZMuYJjNmzfXnXfeOekyJGnd\nePzxx1+oqi1L6bMmA+DOO+9kdnZ20mVI0rqR5G+X2schIEnqKQNAknrKAJCknjIAJKmnDABJ6qlO\nAZBkb5KLSeaSHHmNdv88yY0kv7LUvpKk1bVoACSZAo4B+4C7gPckuWuBdg8Df7nUvpKk1dflCGA3\nMFdVl6rqJeAUsH9Iu98B/gy4toy+kqRV1iUAtgKXW/NXmmWvSLIVeDfwiaX2lSRNxqjuBP4vwL+v\nqh8lWdYLJDkIHATYsWPHiMqStJYdPXf0/0/fd3TBdhqPLgFwFdjemt/WLGubBk41O//NwENJbnTs\nC0BVnQBOAExPT1eX4iVJy9clAM4Du5LsZLDzPgD8WrtBVe28OZ3kD4HPVNWfJ7llsb6SpMlYNACq\n6kaSw8BZYAo4WVUXkhxq1h9fat/RlC5JWolO5wCqagaYmbds6I6/qn5rsb6S+qs97q/J8k5gSeqp\nNfn/ACRtLP7VvzZ5BCBJPWUASFJPGQCS1oSj5446VLTKDABJ6ikDQJJ6ygCQpJ4yACSppwwASeop\nA0CSesoAkKSeMgAkqacMAEnqKR8GJ2lsvLN3bfMIQJJ6ygCQpJ7qFABJ9ia5mGQuyZEh6/cneTLJ\nE0lmk+xprXsuyVM3142yeEnS8i16DiDJFHAMeAC4ApxPcqaqnmk1+wJwpqoqyVuB/wm8pbX+/qp6\nYYR1S5JWqMtJ4N3AXFVdAkhyCtgPvBIAVfX9VvtNQI2ySEn90T5xfPS+owu208p1GQLaClxuzV9p\nlr1KkncneRb4C+B9rVUFPJLk8SQHV1KsJGl0RnYSuKpOV9VbgHcBH2ut2lNVbwP2AR9M8o5h/ZMc\nbM4fzF6/fn1UZUmSFtAlAK4C21vz25plQ1XVo8Cbk2xu5q82n68BpxkMKQ3rd6KqpqtqesuWLR3L\nlyQtV5cAOA/sSrIzya3AAeBMu0GSn0mSZvoe4HXAi0k2JbmtWb4JeBB4epTfgCRpeRY9CVxVN5Ic\nBs4CU8DJqrqQ5FCz/jjwy8B7k7wM/B3wq80VQbcDp5tsuAX4VFV9bkzfiyRpCTo9CqKqZoCZecuO\nt6YfBh4e0u8ScPcKa5S0jvj4h/XDO4ElqacMAEnqKQNAknrKAJCknjIAJKmnDABJ6ikDQJJ6ygCQ\npJ7yfwJLWrN8NPR4eQQgST1lAEhSTxkAktRTBoAk9ZQBIEk9ZQBIUk8ZAJLUUwaAJPVUpwBIsjfJ\nxSRzSY4MWb8/yZNJnkgym2RP176SpMlYNACSTAHHgH3AXcB7ktw1r9kXgLur6m3A+4BPLqGvJGkC\nuhwB7AbmqupSVb0EnAL2txtU1ferqprZTUB17StJmowuAbAVuNyav9Ise5Uk707yLPAXDI4COveV\nJK2+kT0MrqpOA6eTvAP4GPDOpfRPchA4CLBjx45RlSVpFbQf2qb1o8sRwFVge2t+W7NsqKp6FHhz\nks1L6VtVJ6pquqqmt2zZ0qEsSdJKdAmA88CuJDuT3AocAM60GyT5mSRppu8BXge82KWvJGkyFh0C\nqqobSQ4DZ4Ep4GRVXUhyqFl/HPhl4L1JXgb+DvjV5qTw0L5j+l4kbWD+b4DR63QOoKpmgJl5y463\nph8GHu7aV5I0ed4JLEk9ZQBIUk8ZAJLUUwaAJPWUASBJPWUASFJPGQCS1FMGgCT1lAEgST1lAEhS\nTxkAktRTBoAk9ZQBIEk9ZQBIUk8ZAJLUUyP7n8CS+sf/Bby+eQQgST3VKQCS7E1yMclckiND1v96\nkieTPJXki0nubq17rln+RJLZURYvSVq+RYeAkkwBx4AHgCvA+SRnquqZVrNvAD9fVd9Jsg84Adzb\nWn9/Vb0wwrolSSvU5QhgNzBXVZeq6iXgFLC/3aCqvlhV32lmHwO2jbZMSdKodQmArcDl1vyVZtlC\n3g98tjVfwCNJHk9ycKFOSQ4mmU0ye/369Q5lSZJWYqRXASW5n0EA7Gkt3lNVV5O8Efh8kmer6tH5\nfavqBIOhI6anp2uUdUnaWG5efXT0vqMTrWO963IEcBXY3prf1ix7lSRvBT4J7K+qF28ur6qrzedr\nwGkGQ0qSpAnrEgDngV1Jdia5FTgAnGk3SLID+DTwG1X11dbyTUluuzkNPAg8PariJUnLt+gQUFXd\nSHIYOAtMASer6kKSQ83648BHgTcAH08CcKOqpoHbgdPNsluAT1XV58bynUiSlqTTOYCqmgFm5i07\n3pr+APCBIf0uAXfPXy5JmjzvBJaknjIAJKmnDABJ6ikDQJJ6ygCQpJ4yACSppwwASeopA0CSesoA\nkKSeMgAkqacMAEnqKQNAknrKAJCknhrpfwSTtPHd/G9cWv88ApCknjIAJKmnDABJ6qlOAZBkb5KL\nSeaSHBmy/teTPJnkqSRfTHJ3176SpMlY9CRwkingGPAAcAU4n+RMVT3TavYN4Oer6jtJ9gEngHs7\n9pWkZWmfkD5639EF22m4LkcAu4G5qrpUVS8Bp4D97QZV9cWq+k4z+xiwrWtfSdJkdAmArcDl1vyV\nZtlC3g98dql9kxxMMptk9vr16x3KkiStxEjvA0hyP4MA2LPUvlV1gsHQEdPT0zXKuqTXstB17Q4p\naKPrEgBXge2t+W3NsldJ8lbgk8C+qnpxKX2l1eANTNKrdRkCOg/sSrIzya3AAeBMu0GSHcCngd+o\nqq8upa8kaTIWPQKoqhtJDgNngSngZFVdSHKoWX8c+CjwBuDjSQBuVNX0Qn3H9L1Ikpag0zmAqpoB\nZuYtO96a/gDwga59pXFyqEfqxofBSQvw5LA2Oh8FIUk9ZQBIUk8ZAJLUUwaAJPWUJ4G1bnm1j7Qy\nBoC0RH28Osiw3ZgcApKknjIAJKmnDABJ6ikDQJJ6ygCQpJ4yACSppwwASeop7wPQmuc16Oqi/XOy\nke/JGCUDQBqRYUHljkhrmUNAktRTnQIgyd4kF5PMJTkyZP1bkvx1kh8m+fC8dc8leSrJE0lmR1W4\nJGllFh0CSjIFHAMeAK4A55OcqapnWs2+Dfwu8K4FXub+qnphpcVKkkanyxHAbmCuqi5V1UvAKWB/\nu0FVXauq88DLY6hRkjQGXQJgK3C5NX+lWdZVAY8keTzJwYUaJTmYZDbJ7PXr15fw8pKk5ViNk8B7\nquptwD7gg0neMaxRVZ2oqumqmt6yZcsqlCVJ/dYlAK4C21vz25plnVTV1ebzNeA0gyElSdKEdQmA\n88CuJDuT3AocAM50efEkm5LcdnMaeBB4ernFSpJGZ9GrgKrqRpLDwFlgCjhZVReSHGrWH0/yJmAW\n+EngR0k+BNwFbAZOJ7n5tT5VVZ8bz7ei9c47fqXV1elO4KqaAWbmLTvemv4mg6Gh+b4H3L2SAqX1\nbD3/+0gDeePzTmBJ6ikDQJJ6ygCQtOEcPXfUIawODABJ6ikDQJJ6ygCQpJ4yACSppwwASeopA0CS\nesoAkKSeMgAkqac6PQtIGrW+36Sznp8RpI3DIwBJ6ikDQJJ6ygCQpJ4yACSppzoFQJK9SS4mmUty\nZMj6tyT56yQ/TPLhpfSVJE3GolcBJZkCjgEPAFeA80nOVNUzrWbfBn4XeNcy+kpaQ/p+hVafdDkC\n2A3MVdWlqnoJOAXsbzeoqmtVdR54eal9JUmT0SUAtgKXW/NXmmVdrKSvJGmM1syNYEkOAgcBduzY\nMeFqJG0E7eEsb7L7cV2OAK4C21vz25plXXTuW1Unqmq6qqa3bNnS8eUlScvV5QjgPLAryU4GO+8D\nwK91fP2V9NUG4AnFpfEREVpNiwZAVd1Ichg4C0wBJ6vqQpJDzfrjSd4EzAI/CfwoyYeAu6rqe8P6\njuubkSR11+kcQFXNADPzlh1vTX+TwfBOp76SpMnzTmBJ6ikDQJJ6ygCQpJ5aM/cBSJocr9bqJ48A\nJKmnDABJ6ikDQJJ6ygCQpJ7yJLBGwpOI4+UjIjQOHgFIUk95BCCpF3w09I/zCECSesoAkKSeMgAk\nqacMAEnqKU8CSz3lpbsyAKR1zPsDtBKdhoCS7E1yMclckiND1ifJHzTrn0xyT2vdc0meSvJEktlR\nFi9JWr5FjwCSTAHHgAeAK8D5JGeq6plWs33ArubjXuATzeeb7q+qF0ZWtSbGYQNp4+gyBLQbmKuq\nSwBJTgH7gXYA7Af+uKoKeCzJ65PcUVXPj7xiSVohbwob6DIEtBW43Jq/0izr2qaAR5I8nuTgcguV\nJI3WapwE3lNVV5O8Efh8kmer6tH5jZpwOAiwY8eOVShLkvqtyxHAVWB7a35bs6xTm6q6+fkacJrB\nkNKPqaoTVTVdVdNbtmzpVr0kadm6BMB5YFeSnUluBQ4AZ+a1OQO8t7ka6O3Ad6vq+SSbktwGkGQT\n8CDw9AjrlyQt06JDQFV1I8lh4CwwBZysqgtJDjXrjwMzwEPAHPAD4Leb7rcDp5Pc/FqfqqrPjfy7\nkCQtWadzAFU1w2An3152vDVdwAeH9LsE3L3CGjUBXu65vnmDmLrwTmCpRwx2tfkwOEnqKQNAUq8d\nPXe0t0dGBoAk9ZQBIEk95UlgqUfOPXdu6PL77rxvVevQ2mAA9Fxfxz4lGQBSLyz0l7/6zQCQJPr5\niGhPAktST3kE0BOO9ffPUoZ9hrX1xPDG5xGAJPWUASBJPeUQkCTN05cTwgbABuNYv0bFm8Y2PgNA\n2kC83n/0NvLRgAGwjvnXvibBI4ONo9NJ4CR7k1xMMpfkyJD1SfIHzfonk9zTta+klTn33LlXPjRe\nG+3R0YseASSZAo4BDwBXgPNJzlTVM61m+4Bdzce9wCeAezv21SI20g+cRmet7fA9Mlh/ugwB7Qbm\nmv/vS5JTwH6gvRPfD/xx87+BH0vy+iR3AHd26KuGO3ptRBsxGDbKeYEuAbAVuNyav8Lgr/zF2mzt\n2Hfdc8etcVhrf+GP2lK/v7UaGAv9/q+HYFgzJ4GTHAQONrPfT3JxGS+zGXhhdFWN1FquDdZ2fWu5\nNljb9W2Y2v6KvxpjKUOtaNv9Pr8/wlJ+zLDafnqpL9IlAK4C21vz25plXdr8RIe+AFTVCeBEh3oW\nlGS2qqZX8hrjspZrg7Vd31quDdZ2fda2fGu5vlHV1uUqoPPAriQ7k9wKHADOzGtzBnhvczXQ24Hv\nVtXzHftKkiZg0SOAqrqR5DBwFpgCTlbVhSSHmvXHgRngIWAO+AHw26/VdyzfiSRpSTqdA6iqGQY7\n+fay463pAj7Yte8YrWgIaczWcm2wtutby7XB2q7P2pZvLdc3ktoy2HdLkvrGx0FLUk+tuwBI8q+T\nXEjyoyTT89b9XvPIiYtJ/uUC/X8qyeeTfK35/A/HVOf/SPJE8/FckicWaPdckqeadrPjqGWBr3s0\nydVWjQ8t0G7VH+WR5D8mebZ5rMjpJK9foN2qbbuVPA5l3JJsT/K/kzzT/G782yFt7kvy3db7/dFV\nrO8136cJb7ufbW2TJ5J8L8mH5rVZtW2X5GSSa0mebi3rtM9a1u9qVa2rD+AfAz8LnAOmW8vvAr4M\nvA7YCXwdmBrS/z8AR5rpI8DDq1DzfwI+usC654DNE9iOR4EPL9JmqtmObwZubbbvXatQ24PALc30\nwwu9R6u17bpsBwYXQXwWCPB24G9W8b28A7inmb4N+OqQ+u4DPrPaP2dd3qdJbrsh7/M3gZ+e1LYD\n3gHcAzzdWrboPmu5v6vr7gigqr5SVcNuEtsPnKqqH1bVNxhckbR7gXZ/1Ez/EfCu8VQ6kCTAvwH+\ndJxfZ0xeeQxIVb0E3HyUx1hV1V9W1Y1m9jEG949MUpft8MrjUKrqMeDm41DGrqqer6ovNdP/F/gK\ng7vw14uJbbt5fgH4elX97QS+NgBV9Sjw7XmLu+yzlvW7uu4C4DUs9DiK+W6vwT0KMEj728dc178A\nvlVVX1tgfQGPJHm8uRt6Nf1Oc8h9coHDyq7bdJzex+Cvw2FWa9t12Q5rYVuR5E7gnwJ/M2T1zzXv\n92eT/JNVLGux92lNbDsG9ykt9IfapLYddNtnLWsbrplHQbQleQR405BVH6mq/zWqr1NVlWTZl0F1\nrPM9vPZf/3uq6mqSNwKfT/Js81fAir1WfQye2PoxBr+cH2MwTPW+UXzdldZ2c9sl+QhwA/iTBV5m\nbNtuPUry94E/Az5UVd+bt/pLwI6q+n5zvufPGTy9dzWs+fcpgxtVfxH4vSGrJ7ntXmWl+6z51mQA\nVNU7l9GtyyMrAL6V5I6qer45zLy2nBph8TqT3AL8EvDPXuM1rjafryU5zeBQbiS/HF23Y5L/Cnxm\nyKqu23TJOmy73wL+FfAL1QxyDnmNsW27eVbyOJRVkeQnGOz8/6SqPj1/fTsQqmomyceTbK6qsT8n\nqMP7NNFt19gHfKmqvjV/xSS3XaPLPmtZ23AjDQGdAQ4keV2SnQwS+v8s0O43m+nfBEZ2RDHEO4Fn\nq+rKsJVJNiW57eY0g5OfTw9rO2rzxljfvcDXncijPJLsBf4d8ItV9YMF2qzmtlvJ41DGrjnP9N+A\nr1TVf16gzZuadiTZzeB3/8VVqK3L+zSxbdey4JH6pLZdS5d91vJ+V1fjzPYoPxjsrK4APwS+BZxt\nrfsIgzPhF4F9reWfpLliCHgD8AXga8AjwE+NsdY/BA7NW/aPgJlm+s0MztZ/GbjAYPhjtbbjfwee\nAp5sflDumF9fM/8Qg6tKvr5a9TE4gX8ZeKL5OD7pbTdsOwCHbr6/DK5gOdasf4rWFWqrsL32MBjK\ne7K1zR6aV9/hZjt9mcGJ9Z9bpdqGvk9rZds1X38Tgx36P2gtm8i2YxBCzwMvN/u59y+0zxrF76p3\nAktST22kISBJ0hIYAJLUUwaAJPWUASBJPWUASFJPGQCS1FMGgCT1lAEgST31/wCE4jpRMRFVcwAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb7b54210f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "num_bins=50\n",
    "plt.hist(x_res,num_bins,normed=1,facecolor='green',alpha=0.5)\n",
    "plt.hist(y_res,num_bins,normed=1,facecolor='green',alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<bound method MultivariateGaussian._pdf of MultivariateGaussian(\n",
      "    mu=[ 5 -1]\n",
      "    cov=[[ 1.   0.5]\n",
      " [ 0.5  2. ]]\n",
      ")>\n"
     ]
    }
   ],
   "source": [
    "print(rv._pdf())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
